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ABSTRACT 


This  thesis  addresses  optimal  methods  for  the  detection  of  acoustic  signals  corrupted 
by  colored  noise.  In  achieving  this  we  provide  a  study  of  the  characteristics  of  ambient 
noise  in  the  ocean  and  the  digital  techniques  which  can  be  used  in  the  process  of  detecting 
known  acoustic  signals  which  are  corrupted  by  that  noise.  Various  techniques  are  studied, 
in  particular  the  use  of  matrix  decomposition  techniques  applied  to  the  correlation  matrix 
or  to  a  data  matrix,  and  the  matched  filter  for  colored  noise.  Other  methods  such  as  the 
inverse  filter,  the  differential  operator,  and  the  adaptive  prediction-error  filter  will  also  be 
looked  at  for  their  whitening  properties.  The  theoretical  foundations  of  those  techniques  are 
presented  as  well  as  the  application  of  each  method  to  the  problem.  Simulations  are 
conducted  for  each  technique  in  order  to  provide  quantified  performance  measurements 
supporting  the  use  of  each  method. 
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I.  INTRODUCTION 


A.  THESIS  GOAL 

This  thesis  is  concerned  with  the  detection  of  periodic  signals  as  generally  emitted  by 
the  machinery  of  travelling  vessels  whether  submerged  or  at  the  surface  of  the  ocean.  The 
goal  is  to  provide  a  study  of  the  discrete  time  techniques  that  can  be  used  for  the  optimal 
detection  of  acoustic  signals  in  colored  noise.  To  achieve  this  goal,  a  comprehensive  study 
of  the  ambient  noise  present  in  the  ocean  is  discussed  as  well  as  the  effect  on  this  noise  of 
the  recording,  sampling  and  digitization  of  the  input  data.  Having  identified  the 
characteristics  of  the  ambient  noise  presented  to  the  detector  algorithm,  this  thesis  then 
studies  various  transformation  techniques  to  process  the  correlated  samples  received  in 
order  to  achieve  an  optimal  detection  of  the  desired  signal. 

B.  ENVIRONMENT 

The  environment  in  which  underwater  signals  are  generated,  propagated  and 
intercepted  is  of  great  importance  to  the  techniques  of  signal  processing  used  in  analyzing 
propagating  sounds.  We  are  concerned  with  the  underwater  properties  of  the  world’s 
oceans.  The  environment  of  our  oceans  in  terms  of  ambient  noise,  is  characterized  by  the 
propagating  medium,  namely  sea  water.  Water  in  itself  profoundly  changes  the  spectral 
characteristics  of  any  broadband  sound  propagating  in  it,  in  that  it  absorbs  acoustic  energy 
in  a  way  that  is  highly  dependent  on  the  frequency  content  of  those  signals.  Sea  water  is 
particular  in  that  it  contains  significant  amounts  of  dissolved  salts.  Some  of  these  salts  also 
have  important  absorptive  qualities  that  affect  signals  of  various  frequencies  differently. 
These  properties  have  a  significant  effect  on  the  ambient  noise  to  which  all  signals  of 
interest  are  subjected  to.  This  noise  is  quite  different  from  the  common  assumptions  of 
white  noise.  That  is,  the  spectral  shape  of  the  ambient  noise  in  the  ocean  is  frequency 
dependent. 
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C.  DISCRETE  TIME  DETECTION 


The  detection  aspect  of  this  work  is  accomplished  through  the  use  of  probability  and 
statistical  tools  to  design  a  receiver  which  provides  the  user  the  ability  to  discriminate 
between  signal  plus  noise  and  noise  only.  Through  the  use  of  modem  technology,  these 
signals  are  most  often  sampled  in  order  to  be  processed  with  the  digital  processors  in  use 
today.  Therefore,  this  thesis  concentrates  on  discrete  time  detection,  that  is,  detection  using 
only  discrete  time  quantities  as  input  and  ouq)ut  data.  The  theory  of  detection  in  continuous 
time  is  only  addressed  in  those  areas  where  it  is  useful  in  order  to  provide  the  reader  with 
a  greater  understanding  of  discrete  time  detection.  The  receivers  to  be  designed  are 
optimum  in  that  they  best  satisfy  a  given  criteria  such  as  signal  to  noise  ratio,  under  a  given 
set  of  assumptions. 

D.  TECHNIQUES 

Because  of  the  spectral  shape  of  the  ambient  noise  corrupting  the  signals,  the 
techniques  that  are  investigated  are  concerned  with  whitening  the  noise.  This  may  be  done 
in  different  ways.  The  matched  filter  for  colored  noise  inherently  whitens  the  noise  samples 
and  achieves  detection  in  one  step.  Other  optimal  techniques  require  the  use  of  a 
prewhitener  and  a  white  noise  detector.  Since  we  know  the  power  spectral  density 
(spectrum)  of  the  noise,  it  is  possible  to  create  a  filter  which  acts  as  an  equalizer.  Running 
the  received  signal  through  it  has  the  effect  of  whitening  the  additive  noise  allowing  the  use 
of  an  optimal  white  noise  detection  scheme.  This  process  is  known  as  inverse  filu  ring. 
Other  techniques  use  a  decomposed  correlation  matrix  or  covariance  matrix  to  transform 
the  input  data  vector  from  one  consisting  of  both  correlated  signal  and  noise  samples,  to 
one  with  uncorrelated  noise  samples  plus  the  signal.  Various  techniques  have  different 
advantages  in  terms  of  performance,  complexity  and  computational  effectiveness.  The  use 
of  a  data  matrix  to  compute  the  orthogonal  matrices  is  also  addressed  by  using  the  Singular 
Value  Decomposition  (SVD)  and  QR  methods. 


2 


E.  THESIS  OUTLINE 


This  thesis  is  divided  into  two  main  parts.  The  first  part,  addressed  in  Chapter  II, 
provides  a  study  of  the  processes  that  affect  the  detection  of  signals  in  the  underwater 
acoustical  environments.  Acoustic  signals  in  the  ocean  are  necessarily  corrupted  by  the 
ambient  noise  of  the  environment  from  which  they  have  been  recorded  and  in  which  they 
have  been  emitted.  Furthermore,  the  means  in  which  they  are  recorded  further  modifies  the 
inherent  characteristics  of  the  noise  corrupting  the  signal.  In  order  to  provide  an  optimal 
detection  scheme  for  these  signals,  we  must  have  a  good  understanding  of  the  noise 
process,  and  the  transmission  of  the  noise  and  of  the  acoustic  signals  within  the  ocean. 

The  second  major  part  is  dealt  with  in  Chapter  HI.  This  chapter  provides  the  reader 
with  basic  detection  theory  as  applied  to  the  discrete  case  with  colored  noise  as  is  prevalent 
in  the  ocean.  The  major  obstacle  to  be  overcome  for  the  detection  of  signals  in  colored  noise 
is  that  the  samples  of  the  noise  are  correlated.  In  order  for  most  discrete  time  detection 
scheme  to  function,  these  samples  must  undergo  a  transformation  which  decorrelates  the 
noise.  The  various  procedures  in  which  this  decorrelation  can  be  done  will  be  documented. 
We  will  study  the  matched  filter  for  colored  noise  and  prewhitening  of  the  noise  using 
various  factorization  techniques.  The  correlation  matrix,  the  covariance  matrix  or  a  data 
matrix  will  be  factored  into  products  of  other  matrices  which  have  the  required 
characteristics  leading  to  uncorrelated  noise  samples.  We  will  also  address  the  use  of  the 
differential  operator  to  whiten  the  noise,  and,  the  whitening  properties  of  the  prediction- 
error  filter. 

Chapter  IV  presents  the  results  of  simulations  undertaken  with  each  of  the  techniques 
presented  in  Chapter  HI  to  provide  comparative  performance  results.  In  each  case,  white 
noise  discrete  time  detection  of  periodic  signals  is  done  through  the  use  of  the  Fast-Fourier 
Transform  (FFT)  to  compute  averaged  periodograms.  First,  the  averaged  periodogram  is 
shown  for  the  original  sequence  with  the  colored  noise.  Then  the  averaged  periodogram  of 
the  whitened  sequence  is  shown  for  the  particular  transformation  under  examination.  This 
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will  allow  for  an  objective  view  of  the  result  and  permit  us  to  make  qualified  statements  on 
the  performance  of  each  technique. 
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II.  ACOUSTIC  AMBIENT  NOISE  SPECTRUM 


The  spectral  characteristics  of  the  noise  affecting  the  detection  performance  of 
modem  detectors  and  receivers  are  shaped  by  a  variety  of  mechanisms.  Figure  2.1  shows  a 
block  diagram  of  the  various  causes  and  transformations  which  have  a  part  in  shaping  the 
noise  corrupting  received  signals  in  acoustic  environments. 


Processed 

Ambient 

Noise 


Figure  2.1:  Generation  and  processing  of  ocean  acoustic  noise 


A.  NOISE  SOURCES 

The  study  of  the  ambient  noise  frequency  spectrum  begins  with  a  consideration  of  the 
individual  somces  of  noise.  These  sources  are  divided  into  the  following  general  groups: 
Thermal  agitation,  hydrodynamic,  man-made,  seismic,  biological,  and  polar.  The  effect  of 
these  sources  cannot  be  completely  studied  without  considering  the  effects  of  attenuation. 
Attenuation  is  defined  as  the  loss  of  acoustic  energy  from  a  sound  beam.  It  is  divided  in  two 
parts:  absorption  mechanisms  which  convert  acoustics  energy  into  thermal  energy  as  a 
result  of  the  compressive/decompressive  effects  of  the  sound  wave,  and  scattering  effects 
which  deflect  energy  out  of  the  direction  of  propagation  of  the  acoustic  wave  [Ref.  1]. 

1.  Thermal  Agitation  Noise 

The  minimum  noise  level  of  a  medium  is  determined  by  the  effects  of  thermal 
agitation.  For  the  ocean,  in  the  temperature  range  of  0-30*  C,  the  equivalent  thermal-noise 

sound  pressure  level  (SPL)  in  dB  re  2  x  lO^pPa  for  a  IHz  bandwidth  is  given  by  [Ref.  2] 

SPLjfj »  -  101  +  201og/ .  (2.1) 

The  thermal  noise  spectral  density  derived  from  this  equation  is  shown  in  Figure  2.2. 
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Figure  2.2:  Thermal  agitation  sound  pressure  level  (SPL) 

According  to  this  equation,  the  thermal  noise  spectrum  has  a  level  of  -10  dB  at  35  kHz 
and  a  positive  slope  of  +6  dB/octave.  At  higher  frequencies,  20-30  kHz  and  up,  the 
observed  minimum  ambient  noise  levels  are  very  close  to  the  thermal  noise  levels  shown 
on  Figure  2.2.  At  lower  frequencies  (/  <  10  kHz)  however,  a  significant  difference  exists 
where  the  observed  ambient  noise  levels  are  significantly  higher  than  those  which  would 
result  only  from  thomal  noise. 

2.  Hydrodynamic  Noise 

This  group  contains  the  noise  sources  which  are  related  to  hydrodynamic  processes. 
Many  of  these  are  important  radiators  of  sound  and  contribute  to  the  noise  spectra  in  the 
ocean. 

a.  Air  Bubbles 

Air  bubbles  affect  the  ambient  noise  spectrum  in  two  distinct  ways.  Their  efiect 
on  the  attenuation  properties  of  the  medium  are  discussed  in  a  later  paragraph.  This 
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paragraph  addresses  the  manner  in  which  soimd  energy  is  radiated  when  air  bubbles  exist, 
and  the  cavitation  noise  resulting  from  their  collapse.  When  wind  effects  are  considered, 
substantial  noise  levels  [Ref.  3]  have  been  observed  when  air  is  entrained  into  the  water  by 
water  droplets  on  the  surface.  However,  bubbles  do  not  only  occur  as  a  result  of  wind. 
Bubbles  are  also  created  as  a  result  of  biological  effects  such  as  decaying  matter,  sea  floor 
gas  seepage,  fish  and  marine  mammals.  As  air  bubbles  in  the  ocean  are  transported  or 
simply  rise  to  the  surface,  they  are  subjected  to  temperature  and  pressure  fluctuations  due 
to  depth,  turbulence,  currents  and  waves.  These  variations  induce  oscillations  in  the 
bubbles  causing  them  to  emit  noise  at  frequencies  related  to  their  diameters.  It  has  been 
concluded  [Ref.  3]  that  air  bubbles  and  their  cavitation  produced  at  or  near  the  surface  as  a 
result  of  the  wind,  are  an  important  source  of  wind  dependent  noise  in  the  frequency  range 
of  50-10,000  Hz.  The  spectrum  shape  of  cavitation  noise  is  similar  to  that  of  the  air  bubble 
noise,  and  is  characterized  by  a  slope  of  approximately  -6  dB/octave  in  the  frequency  range 
identified  above.  As  the  frequency  is  reduced,  the  noise  level  resulting  firom  this  effect 
drops-off  with  a  12  dB/octave  slope  due  the  lack  of  bubbles  of  a  size  large  enough  to 
resonate  at  those  frequencies.  Experimentation  has  shown  that  at  shallow  depth,  both  the 
noise  spectrum  level  and  its  shape  in  fr-equencies  between  50-10,000  Hz  are  likely  to  be 
mostly  the  result  of  resonating  noise  and  cavitation  noise  from  the  air  bubbles  [Ref.  3], 

b.  Surface  Waves 

The  performance  of  an  underwater  transducer  is  affected  by  the  subsurface 
pressure  fluctuations  caused  by  changes  in  surface  elevations.  In  the  frequency  range  of 
500  Hz  to  20  kHz,  the  agitation  of  the  sea  surface  is  the  primary  source  of  ambient  noise 
[Ref.  1].  Relations  between  sea  state,  mean  wave  height,  and  representative  wind  speed 
have  been  derived  and  are  usually  given  in  table  formats  such  as  that  giv«i  at  [Ref.  1:  p. 
414].  As  a  result,  it  is  possible  to  characterize  the  ambient  noise  caused  by  surface  waves, 
by  specifying  wind  speed.  In  the  range  of  frequencies  specified  above,  it  has  been 
established  that  the  noise  spectrum  level  decreases  at  a  rate  of  17  dB/decade  [Ref.  1]. 
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c.  Turbulence 

Turbulence  is  defined  by  the  instability  of  the  water  motion.  It  is  usually  created 
at  boundary  regions  where  water  and  solid  surfaces  interact  or  when  fluid  layers  flow  past 
one  another  with  different  velocities  and  directions.  Turbulence  between  the  transducer  or 
the  vessels  containing  it  is  generally  considered  to  be  self-noise.  Turbulent  effects  usually 
consist  of  large  volume  fluctuations  resulting  in  pressure  variations  at  frequencies  less  than 
1  Hz.  Although  some  noise  is  created  at  higher  frequencies  when  large-scale  eddies 
breakup  into  smaller  scale  phenomena,  it  has  been  observed  [Ref.  3]  that  these  noise  levels 
are  several  orders  of  magnitude  below  the  ambient  noise  level  observed.  It  is  considered 
that  the  radiated  noise  caused  by  turbulence  is  not  a  contributing  factor  to  the  ambient  noise 
level  in  the  ocean.  However,  the  pressure  fluctuations  do  affect  pressure  transducers  when 
located  in  a  region  of  turbulence.  Substantial  ambient  oceanic  turbulence  has  been  observed 
[Ref.  3],  with  pressure  level  spectrum  occurring  in  the  frequency  range  below  10  Hz.  In 
swift  oceanic  streams,  the  pressure  level  can  be  affected  up  to  100  Hz,  while  extreme  tidal 
currents  can  generate  pressures  at  even  higher  frequencies.  It  has  been  concluded  that  the 
turbulent  pressure  fluctuations  are  an  important  component  of  noise  below  10  Hz,  and  in 
some  cases  below  100  Hz,  but  that  its  contribution  to  the  ambient  noise  spectrum  from 
radiated  noise  is  negligible  [Ref.  3]. 

3.  Man-made  Sources 

Ships  are  important  sources  of  underwater  sounds.  At  any  given  time,  it  is  estimated 
that  there  are  over  1000  ships  under  way  in  the  North  Atlantic  Ocean  [Ref.  4].  Rotational 
and  reciprocating  machinery  is  required  for  propulsion,  control,  and  habitability.  These 
components  generate  vibrations  which  are  transmitted  through  the  hull  of  the  ship  to  the 
surrounding  water.  The  propeller  of  the  ships  is  a  major  source  of  this  underwater  noise.  At 
lower  frequencies  (10  Hz  to  250  Hz)  than  those  associated  with  surface  waves,  a  major 
source  of  noise  is  generated  by  distant  shipping.  The  sh^  of  the  noise  spectra  and  its 
strength  resulting  from  this  oceanic  traffic  will  be  affected  by  a  combination  of  the 


8 


transmission  loss,  the  number  of  ships,  and  the  distribution  of  those  ships.  For  example,  a 
transducer  in  deep  water  in  an  open  ocean  may  receive  significant  shipping  noise  from 
widely  scattered  ships  because  the  average  transmission  loss  will  be  relatively  small.  The 
ambient  noise  may  be  likewise  considerably  affected  in  a  situation  where  transmission 
losses  are  high  but  ship  concentration  is  large,  in  a  shallow  water  location  such  as  that 
occurring  along  coastal  shipping  lanes.  Traffic  noise  may  also  depend  on  the  types  of  ships 
involved,  each  of  which  may  generate  different  noises.  However,  we  assume  that  the 
number  of  ships  affecting  the  typical  acoustic  location  will  blend  the  individual  differences 
of  the  distant  ships  into  an  average  source  characteristics  with  Gaussian  distribution. 

For  most  surface  ships,  the  effective  source  of  the  radiated  noise  is  the  engine/ 
propeller  combination  where  the  noise  is  generated  between  three  and  ten  meters  below  the 
surface.  This  places  the  source  in  the  shallow-water  chaimel.  As  a  result  of  the  surface 
reflection,  the  source  and  its  image  will  operate  as  an  acoustic  doublets  radiating  noise  with 
a  spectrum  slope  of  +6  dB/octave  relative  to  the  spectrum  of  the  simple  source. 

Studies  of  noise  from  surface  ships  have  determined  [Ref.  3]  that  the  sound  pressure 
level  spectrum  has  a  slope  of  about  -6  dB/octave.  However,  the  spectrum  is  highly  variable 
below  1  kHz,  and  some  tests  have  shown  a  general  levelling  of  the  spectra  at  the  100  Hz 
frequency. 

Industrial  activities  may  also  locally  affect  the  ambient  noise  characteristics.  Many 
mechanical  movements  related  to  such  things  as  oil  drilling,  pile  driving,  mining,  etc.,  will 
transmit  sounds  to  the  water  when  located  in  an  area  close  to  shore  or  at  sea  such  as  on 
drilling  platforms.  These  noises  can  be  accounted  for  in  the  design  of  fixed  transducers  in 
their  vicinity  but  are  generally  ignored  in  the  design  of  other  transducers  due  to  their  highly 
localized  nature.  Operators  however  must  be  aware  of  these  sources. 

4.  Seismic  Noise 

The  earth’s  crust  is  constantly  moving  as  a  result  of  volcanic  and  tectonic  action.  This 
action  produces  waves  which  are  transmitted  through  the  solid  crust  and  eventually  to  the 
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water.  Due  to  the  nature  of  the  crust  and  of  the  water,  even  far  from  the  point  of  origin  of  a 
disturbance,  appreciable  amounts  of  energy  can  find  their  way  to  the  ocean  waters  and  be 
propagated  as  compressions!  waves.  The  same  effect  can  results  from  man-made 
explosions.  The  spectral  characteristics  of  such  an  event  wiU  depend  on  the  magnitude  of 
the  disturbance,  the  range,  and  the  propagation  path.  This  noise  is  typically  characterized 
by  one  or  a  series  of  transients  of  relative  short  duration.  In  general,  the  affected  frequencies 
will  be  less  than  500  Hz  with  a  maximum  between  2  and  20  Hz  [Ref.  3].  Noticeable 
underwater  sound  will  be  generated  from  1  to  100  Hz  [Ref.  3].  Outside  of  the  specific  effect 
of  the  transients  themselves,  a  background  of  continuous  seismic  noise  is  also  present.  This 
effect  is  attributed  to  the  after  effects  of  the  main  transients  and  to  the  effects  of  storms,  and 
waves  and  swells  hitting  the  shores  along  coastal  zones.  It  was  found  that  seismic  noise 
dominates  the  ambient  noise  spectrum  at  very  low  frequencies  between  1  and  100  Hz  [Ref. 
3].  However  this  noise  spectrum  level  is  of  a  transitory  nature,  and  is  highly  dependent  on 
time  and  location.  Due  to  its  nature,  seismic  noise  is  particularly  important  to  the  design 
and  operation  of  bottom-mounted  transducers. 

5.  Biological  Sources 

Underwater  sounds  originating  from  marine  life  can  be  an  important  source  of  noise 
at  most  frequencies  of  interest  [Ref.  3].  The  wide  variety  of  sounds  are  usually  of  short 
duration  but  are  often  repeated.  Continuous  noise  is  often  observed  when  it  is  the  result  of 
a  great  number  of  individuals  such  as  shrimps.  As  with  most  other  sources,  the  ambient 
noise  of  biological  origin  varies  with  location,  time  and  species  involved.  In  some  case,  the 
sounds  are  correlated  with  geographical,  seasonal  and  diurnal  pattern  enabling  users  of 
systems  to  improve  the  performance  of  their  systems  by  considaing  the  effect  of  these 
sources.  However,  because  of  this  variance,  it  is  generally  ignored  in  establishing  the 
ambient  noise  spectrum  of  an  area  of  ocean. 
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6.  Polar  Noise 

Polar  regions  cause  noise  spectra  which  are  very  different  than  those  prevalent  in  other 
oceans  because  of  the  effect  of  the  ice  covering  large  portions  of  the  water.  Cracking  of  the 
ice  as  a  result  of  thermal  effects,  and  movements  of  large  pieces  generate  cracking, 
straining,  crunching,  sliding  and  percussion  sounds.  Other  sounds  are  created  by  icebergs 
as  they  melt.  As  the  ice  changes  into  water,  numerous  pressurized  air  bubbles  which  were 
trapped  in  the  ice,  burst  through  the  thinning  surface.  As  the  air  rushes  out,  a  sharp  cracking 
sound  is  produced.  The  noise  produced  by  these  pressurized  air  bubbles  has  been  formd  to 
be  the  prevalent  source  of  noise  when  icebergs  are  present.  These  sounds  are  highly 
characteristic  and  must  be  addressed  differently  than  more  typical  underwater  noises  due 
to  their  impulsive  nature.  Some  methods  dealing  with  this  enviroiunent  have  been  devised 
and  generally  include  a  provision  for  removing  the  impulsive  components  of  the  ambient 
noise  prior  to  the  whitening  and  detection  process  [Ref.  5].  However,  the  level  and 
character  of  the  noise  spectra  vary  according  to  ice  condition,  wind  speed,  ice  pack  snow 
cover,  and  air  temperature  changes  [Ref.  4].  Due  to  the  extreme  variations  observed,  it  is 
not  possible  to  predict  the  spectra  and  amplitude  distribution  of  the  noise  to  be  expected  in 
polar  regions. 
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7.  Estimated  Noise  Spectra 

The  nominal  shape  of  the  noise  spectrum  level  for  the  deep-water  case  is  shown  in 
Figure  2.3. 


Figure  2.3:  Deep-water  ambient  noise  spectrum  level.  From  Ref.  [1]. 


Below  20  Hz,  it  has  been  foimd  that  ocean  turbulence  and  seismic  noise  predominate. 
Above  20  Hz  and  until  about  300  Hz,  shipping  noise  and  biological  noises  are  the  nmifi 
contributor.  From  300  Hz  and  above,  surface  noise  highly  dq^endent  on  wind  speed 
provides  most  of  the  noise  energy.  Although  not  shown  on  this  figure,  at  50  kHz,  thermal 
agitation  of  the  water  molecules  overwhelms  all  other  sources  of  noise  and  continues  to 
increase  with  frequency  at  a  rate  of  6  dB/octave  as  shown  in  Figure  2.2. 
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8.  Distribution  of  the  Noise 

It  has  been  shown  that  the  ambient  noise  spectra  in  the  oceans  is  the  result  of  a  very 
large  number  of  noise  sources  of  various  characteristics  whose  sounds  propagate  in  and  are 
affected  by  the  underwater  medium.  Although  some  impulsive  sources  of  underwater 
noises  can  be  non-Gaussian  and  others  may  even  have  asymmetric  distributions,  by  the 
central  limit  theorem,  it  is  reasonable  to  assume  that  in  general  the  noise  distribution  will 
be  Gaussian  in  nature. 

B.  SOUND  ABSORPTION  IN  SEA  WATER 

Losses  of  energy  in  a  sound  wave  can  be  grouped  into  two  types,  spreading  losses  and 
attenuation  losses.  Spreading  losses  results  from  the  simple  geometric  effect  resulting  from 
waves  emitted  from  sources  projecting  in  all  directions.  These  kinds  of  losses,  although 

great  (spherical  spreading  losses  proportional  to  I /range  ,  cylindrical  spreading  losses 
proportional  to  1/ range )  do  not  affect  the  spectral  characteristics  of  the  ambient  noise 
because  aU  frequencies  are  equally  attenuated  by  spreading  losses.  The  major  loss  affecting 
the  noise  spectrum  is  the  attenuation  loss. 

1.  Absorption  Coefilcient  of  Sea  Water 

Any  medium  in  which  sound  propagates  is  absorptive.  This  means  that  part  of  the 
sound  wave  energy  is  lost  through  heat  as  the  wave  propagates  through  the  medium.  By 
convention,  the  attenuation  loss  is  identified  as  a  and  is  expressed  in  dB  per  kilometer. 
What  makes  a  so  important  in  characterizing  the  ambient  noise  spectrum  is  that  it  has  a 
complicated  dependence  on  frequency.  The  complication  arises  because  different  loss 
processes  dominate  different  portion  of  the  frequency  spectrum.  As  a  result  of  these 
complications  most  of  the  current  knowledge  of  the  attenuation  coefficient  comes  from 
measured  data  taken  from  various  section  of  the  ocean  and  from  laboratory  measurements. 
As  a  result  of  the  complexity  of  the  processes  and  of  the  changing  nature  of  the  oceans,  no 
theoretical  curves  has  been  found  to  define  the  attenuation  losses  in  a  specific  area,  rather. 
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scientists  have  relied  on  curve  fits  over  various  parts  of  the  frequency  spectrum  to  define 
the  magnitude  of  the  attenuation  coefficient.  Three  main  effects  cause  the  absorption  of 
sound  in  sea  water.  The  first  is  the  result  of  shear  viscosity.  This  effect  was  studied  by 
Rayleigh  who  derived  an  expression  equivalent  to 


a  = 


1671 V 

3pc^ 


where  a  is  the  intensity  absorption  coefficient, 


-7^  (2.2) 

is  the  shear  viscosity,  p  is  the  density,  c 


is  the  sound  velocity,  and  /  is  the  frequency  (Hz)  [Ref.  6]. 

The  value  of  the  absorption  coefficient  calculated  using  this  formula  which  considers 
only  shear  viscosity,  is  only  about  one  third  of  the  attenuation  actually  observed  in  pure 
distilled  water.  A  second  kind  of  viscosity  called  volume  viscosity  adds  to  the  value  of  the 
absorption  coefficient.  This  occurs  as  a  result  of  the  time  lag  required  for  the  water 
molecules  to  flow  into  lattice  interstices  in  the  crystal  structure.  Because  of  this  effect,  the 
absoiption  coefficient  formula  becomes 


where  is  the  volume  viscosity  coefficient  [Ref.  4]. 


(2.3) 


Figure  2.4  presents  three  different  curves  for  the  value  of  the  absorption  coefficient  as 
a  function  of  frequency.  The  first  curve  shows  the  measured  values  for  sea  water  talcen  from 
Urick  [Ref.  4].  The  second  curve  shows  the  measured  absorption  coefficient  in  distilled 
water  which  considers  both  the  shear  and  the  volume  viscosity  factors  in  accordance  with 
Equation  (2.3).  The  third  shows  the  theoretical  absorption  using  only  the  shear  viscosity 
calculated  with  Equation  (2.2). 
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Figure  2.4:  Comparison  of  absorption  coefficients  in  sea  water,  due  to  shear  and  volume 

viscosity,  and  shear  viscosity  al(»e. 

It  is  evident  from  the  figure  that  not  all  absorption  is  accounted  for  by  using  only  the 
effects  of  shear  viscosity  and  volume  viscosity.  This  is  particularly  true  in  frequencies 
below  100  kHz  which  are  of  particular  interest  for  the  siMiar  detection  problems.  For 
example,  the  absorption  coefficient  of  sea  water  in  the  frequency  range  of  5  to  50  kHz  is 
actually  about  30  times  that  of  distilled  water  [Ref.  4].  It  has  been  found  that  in  these 
frequencies  the  additional  absorption  mechanisms  are  due  to  the  diffo-ent  salts  dissolved  in 
sea  water.  Leonard,  Combs,  and,  Skidmore  [Ref.  7]  have  shown  that  in  addition  to  water, 
two  other  components  of  sea  water  have  been  found  to  be  inqxxiant  to  the  strength  of  the 
absoiption  coefficient.  These  are  magnesium  sulfate  (MgS04)  and  boric  acid  (B(OH)3).  In 
particular,  in  the  frequencies  of  most  interest  for  the  detectimi  problem,  10  to  5  kHz,  the 
effect  of  both  absorptive  components  is  very  important.  TTie  principal  component  of  the 
salts  dissolved  in  sea  water.  Sodium  Chloride  (NaCl),  does  not  appear  to  be  a  significant 
absorptive  constituent.  Despite  the  fact  that  MgS04  constitutes  only  about  4.7  percent  by 
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weight  of  the  total  dissolved  salts  in  sea  water,  it  has  been  found  to  be  the  dominant 
absorptive  constituent  in  sea  water  for  frequencies  below  100  kHz  [Ref.  7].  The  absorption 
is  due  to  the  ionic  relaxation  of  the  MgS04  molecules  in  sea  water.  This  is  a  process  during 
which  MgS04  ions  dissociate  and  reassociate  during  a  finite  time  interval,  called  the 
relaxation  time,  due  to  the  pressure  effect  of  a  sound  wave.  Therefore,  it  is  considered  that 
soimd  absorption  in  sea  water  is  the  sum  of  the  absorption  due  to  water,  magnesium  sulfate, 
and  boric  acid. 

Fisher  and  Simmons,  in  [Ref.  8]  show  that  the  absorptive  coefficient  also  has  a 
dependence  on  pressure.  However  since  pressure  does  not  vary  with  frequency,  but  rather 
with  depth,  its  effect  on  the  spectral  shape  of  the  ambient  noise  is  not  relevant  to  this  study. 
Figure  2.5  shows  the  effect  of  magnesium  sulfate  and  boric  acid  on  the  absorption 
coefficient  of  fresh  water.  These  curves  were  derived  from  laboratory  measurements  [Ref. 
8]. 


Figure  2.5:  Sound  absorption  at  5°C,  1  atm,  35  ppt  salinity,  and  pH=8.0.  After  Fisher  and 

Simmons  [Ref.  8]. 
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2.  Resonating  Bubbles 

Gas  bubbles  in  the  oceans  can  generate  extremely  high  sound  attenuation.  Although 
the  concentration  of  bubbles  in  the  main  body  of  the  ocean  is  insignificant  compared  to  the 
other  sources  of  attenuation,  breaking  waves,  schools  of  fishes  with  their  internal  air  sacs, 
and,  travelling  ships  and  submarines  often  create  large  concentrations  of  gas  bubbles.  A 
sound  wave  passing  through  an  area  with  many  gas  bubbles  will  lose  energy  as  a  result  of 
the  viscous  forces  and  heat  conduction  generated  by  the  compressive/decompressive 
effects  of  the  wave.  The  presence  of  the  gas  bubbles  affects  the  physical  characteristics  of 
the  medium  by  changing  its  sound  speed  and  density.  This  will  cause  reflective  and 
refractive  effects  which  will  further  dissipated  the  energy  of  the  sound.  Finally,  since  the 
gas  bubbles  in  the  water  create  an  inhomogeneous  medium,  the  sound  wave  will  also  be 
scattered,  meaning  that  some  energy  will  be  redirected  in  all  directions  as  a  result  of 
encountering  a  bubble.  The  scattering  effects  is  largely  dependent  on  the  restmant 
frequency  of  the  bubbles.  The  resonant  frequency  is  given  by 


(2-4) 


where  a  is  the  bubble  radius,  y  is  the  ratio  of  specific  heat  for  the  gas  in  the  bubble,  is 
the  hydrostatic  pressure  at  the  depth  of  the  bubble,  and,  is  the  density  of  the  surround¬ 
ing  medium  [Ref.  1]. 

Frequencies  in  the  acoustic  wave  at  or  close  to  the  resonant  frequencies  of  the  bubbles 
can  have  great  influence  on  the  scattering  effect  while  those  frequencies  which  are  fuilher 
from  the  resonant  frequency  may  not  be  scattered.  Fish  detectors  can  find  schools  offish 
and  determine  the  size  of  the  fishes  present  by  extrapolating  their  size  from  the  HnmiTiiwif 
resonant  frequency  o-eated  by  their  air  sacs.  The  measure  of  attenuation  of  a  bubble  is 
determined  by  its  effective  cross  section  a .  This  value  is  a  measure  of  the  fraction  of  enorgy 
that  the  bubble  will  extract  from  a  sound  beam  of  1  square  meter  cross  section,  ff  an 
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acoustic  wave  contains  only  frequencies  at  the  resonant  frequency  f  ^  of  the  bubble,  the 

effective  cross  section  can  be  more  than  a  thousand  times  greater  that  its  actual  cross 
section  [Ref.  1].  When  the  frequencies  of  the  acoustic  wave  are  greater  than  the  resonant 
frequency  of  the  bubble,  the  effective  cross  section  of  the  bubble  affecting  the  wave  is 
approximately  the  same  as  that  of  the  bubbles.  When  the  frequencies  are  smallw,  the 
effective  cross  section  is  much  less  than  the  actual  cross  section  of  the  bubble  and  quickly 
becomes  negligible.  Since  bubbles  rarely  exist  all  with  the  same  radius,  those  bubbles  with 
the  greatest  scattering  effect  will  be  those  whose  resonant  frequencies  are  closest  to  the 
frequencies  contained  in  the  acoustic  wave.  Since  ambient  noise  contains  all  frequencies, 
the  presence  of  bubbles  in  the  area  of  interest  will  further  change  the  frequency  spectra  of 
the  ambient  noise.  However,  there  is  a  practical  limit  to  the  size  of  the  bubbles  in  the  ocean. 
Ocean  conditions  resulting  in  the  creation  of  bubbles  can  be  expected  to  generate  bubbles 
with  a  distribution  around  one  predominant  bubble  size.  It  can  also  be  expected  that  bubbles 
can  only  attain  a  certain  maximum  size  and  therefore  will  only  affect  frequencies  which  are 
above  the  resonant  frequency  of  the  largest  bubbles.  Therefore,  the  ambient  noise  spectrum 
in  regions  affected  by  bubbles  will  have  a  maximum  at  the  frequency  associated  with  the 
predominant  bubble  size.  Its  shape  will  depend  on  the  distribution  of  the  bubbles  size  and 
the  amplitude  of  oscillations  in  the  affected  frequencies. 

C.  TRANSDUCER  EFFECTS 

The  particular  environment  in  which  a  receiver  is  to  used  is  of  prime  importance  to  its 
designers.  For  the  design  of  transducers,  the  most  important  considerations  arise  from  the 
physical  properties  of  water.  In  addition  to  these  considerations,  transducers  must  be 
designed  to  withstand  the  effects  of  sea  water,  biological  activity,  and,  often  large 
hydrostatic  pressures  [Ref.  9]. 

The  ultimate  performance  of  a  transducer  is  limited  by  the  signal  to  noise  ratio  of  its 
output.  The  noise  level  at  this  point  is  the  result  of  the  ambient  noise  in  the  ocean  at  die 
receiving  location  and  of  the  additional  noise  caused  by  the  transducer  itself.  Turbulent 


18 


boundary  layer  pressure  fluctuations,  mechanical  vibrations,  and,  electrical  noise,  are  but 
some  of  the  possible  causes  of  noise.  Obviously  the  success  of  the  design  will  be  in  large 
part  dependent  on  how  well  these  noise  sources  have  been  dealt  with.  Each  of  these  causes 
may  affect  different  portions  of  the  frequency  spectrum  based  on  their  inherent 
characteristics.  Hydrophones  also  have  an  inherent  sensitivity  which  vary  with  frequency. 
Many  different  materials  are  used  in  the  design  of  modem  transducers;  quartz,  piezoelectric 
polymers  such  as  polyvinylidene  (PVDF),  composite  ceramics,  etc.  Each  of  these  materials 
exhibit  different  characteristics  and  will  have  different  frequency  responses.  In  the  design 
of  the  transducer,  these  characteristics  are  taken  into  consideration  and  the  material  selected 
is  usually  used  in  a  spectral  region  where  its  frequency  response  is  relatively  constant  and 
its  sensitivity  adequate  for  its  mission.  When  the  frequency  response  is  not  constant,  it  must 
then  be  “equalized”  by  means  of  calibration  curves  derived  experimentally  by  the  user. 
Nevertheless,  this  experimental  “equalization”  necessarily  results  in  some  distortion  of  the 
received  signal.  Notwithstanding  the  material  chosen,  many  transducers  perform  very 
poorly  at  very  low  frequency  (<  10  Hz)  often  as  a  result  of  having  to  be  shielded  from  large 
hydrostatic  pressures. 

D.  DISCRETE  TIME  SAMPLING 

The  sampling  process  is  the  operation  by  which  an  analog  signal  is  converted  into  a 
corresponding  sequence  of  discrete  time  samples  which  are  uniformly  spaced  in  time. 
Tranducers  are  inherently  analog  sensors.  The  pressure  information  they  convey  must  be 
processed  by  an  acoustic  processor.  The  modern  acoustic  processor  is  invariably  a  digital 
computation  device,  therefore  requiring  discretization  of  the  captured  imderwater  signals. 
The  derivation  of  the  sampling  theorem  is  presented  in  most  books  on  electrical 
engineering  and  will  not  be  repeated  here.  The  sampling  theorem  for  band-limited  signals 
of  finite  energy  can  be  stated  as  follows  [Ref.  8]: 
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“A  band-limited  signal  of  finite  energy,  which  has  no  components  higher  than  W  Hz, 
is  completely  described  by  specifying  the  values  of  the  signal  at  instants  of  time  separated 
by  less  then  1I(2W)  seconds.” 

The  sampling  rate  of  2W  samples  per  second  is  known  as  the  Nyquist  rate.  The 
corollary  of  this  theorem,  is  that  if  the  signal  to  be  sampled  contains  frequency  components 
higher  than  W  Hz  then  it  will  not  be  completely  described  by  its  sampled  values.  This  will 
result  in  a  distortion  of  the  signal  known  as  aliasing.  This  factor  refers  to  the  effect  wherein 
components  of  the  signal  which  are  at  a  frequency  above  W,  wUl  instead  appears  at  a 
different  frequency  in  the  spectrum  of  the  sampled  signal.  The  distortion  caused  by  aliasing 
can.  be  dealt  with  by  two  methods.  Aliasing  can  be  eliminated  by  sampling  the  sigr  .al  at  a 
greater  frequency  than  the  K  yquist  rate  of  the  highest  frequency  component  present  in  the 
signal.  However  when  frequency  components  are  present  with  frequencies  which  are  above 
the  capabilities  of  the  sampling  system,  an  analog  low-pass  filter  can  be  appUed  to  the 
signal  prior  to  sampling,  to  remove  all  frequencies  above  one  half  of  the  sampling  rate.  This 
last  method  is  usually  favored  due  the  availability  of  good,  inexpensive  low-pass  filters. 

Nevertheless,  this  process  necessarily  changes  the  power  spectrum  of  the  noise.  For 
example  when  white  noise  is  corrupting  the  signal.  This  noise  once  filtered  no  longer  has 
equal  power  for  all  frequencies,  but  rather  becomes  a  band-pass  white  noise  process  which 
has  constant  power  within  a  spectral  region  less  than  one  half  the  sampling  frequency  and 
is  zero  outside.  Since  perfect  filters  do  not  exist  in  nature,  the  filter  will  also  exhibit  some 
distortion  especially  near  its  cut-off  frequency  which  will  affect  the  spectral  representation 
of  the  noise  being  passed  through  it. 

E.  NOISE  POWER  SPECTRAL  DENSITY 

The  ambient  noise  power  spectral  density  is  the  result  of  each  of  the  effects  discussed 
above.  Of  course  the  main  contribution  to  the  shape  of  the  ambient  noise  spectra  is  the 
attenuation  of  the  sound  waves  by  sea  water.  Nevertheless,  the  noise  sources,  their  number, 
intensity,  distribution  and  type,  also  have  an  important  effect  on  the  form  of  the  spectral 
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shape  of  the  ambient  noise.  Figure  2.6  summarizes  the  characteristics  of  the  ambient  noise 
spectrum  (i.e.,  the  terms  spectrum  and  spectral  density  are  used  interchangeably)  resulting 
from  the  noise  sources  and  mechanism  discussed  in  this  chapter. 


Figure  2.6:  Ambient  noise  spectrum  composite.  From  Ref.  [3]. 

Figure  2.7  shows  the  ambient  noise  spectra  of  the  ocean  from  the  Point  Sur  SOSUS 
array.  The  signal  was  sampled  at  a  frequency  of  8  kHz  resulting  in  a  mflyiTnnm  frequency 
of  4  kHz. 
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Figure  2.7:  Normalized  average  sound  pressure  level  of  ambient  noise  at  Point  Sur  SOSUS 

array  (sampling  frequency  of  8  kHz) 

As  can  be  observed,  the  noise  sound  pressure  level  of  the  Point  Sur  array  conforms  to 
the  general  shape  shown  in  Figure  2.6.  At  very  low  frequencies  (10-100  Hz),  shipping  noise 
is  clearly  shown.  At  a  frequency  of  approximately  200  Hz  the  characteristics  hump  caused 
by  the  addition  of  the  wind  dependent  noise  to  the  remaining  high  frequency  components 
of  the  shipping  noise  is  observed.  At  the  higher  frequencies  shown  on  this  figure,  surface 
waves  and  agitation  resulting  from  wind  effects  forms  the  general  shape  of  the  curve.  In 
this  case,  the  wind  dependent  noise  appears  to  be  low  and  the  slope  of  the  curve 
corresponds  approximately  to  the  lower  limit  of  prevailing  noise  shown  on  Figure  2.6. 
From  Figure  2.7,  it  can  also  be  seen  that  in  the  frequency  range  shown,  the  slope  drops-off 
approximately  52  dB  for  two  decades  or  at  a  rate  26  dB/decade. 
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III.  OPTIMAL  DETECTION  IN  COLORED  NOISE 

The  concept  of  optimal  detection  addresses  the  problem  of  conceiving  the  best 
detector  at  meeting  one  or  several  criteria  given  a  number  of  assumptions.  An  example  of 
a  common  criteria  used  for  detection  is  the  signal  to  noise  ratio  (SNR).  Therefore,  to  be 
considered  optimal  in  the  sense  of  the  SNR,  a  detector  must  provide  the  maximum  SNR 
possible  given  the  assumptions. 

There  are  two  general  approaches  in  deriving  an  optimal  detector  for  known  signals  in 
colored  noise.  Optimal  solutions  for  signals  in  white  noise  are  well  known  and  attractive  in 
their  simplicity.  Therefore,  if  a  linear  transformation  can  be  applied  to  the  received  signal 
which  whitens  the  noise  while  leaving  the  desired  signal  component  intact,  then  this 
prefiltering  operation  can  be  combined  with  an  optimal  white  noise  detector  as  shown  in 
Figure  3.1. 
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T 
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Figure  3.1:  Optimal  detection  using  whitening  filter 


Several  methods  are  addressed  in  this  chapter  which  can  be  used  to  whiten  the  colored 
noise.  The  noise  is  assumed  to  be  Gaussian. 

Alternatively,  an  optimal  detector  can  also  be  derived  directly  in  the  form  of  a  matched 
filter  for  colored  noise.  This  filter  has  the  advantage  of  removing  the  prewhitening  step  in 
the  detection  process  and  is  shown  in  Figure  3.2. 
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Figure  3.2:  Optimal  colored  noise  detector 
A.  DISCRETE  TIME  DETECTION 

Since  the  first  form  of  optimal  detector  for  colored  noise  uses  a  white  noise  detector, 
it  is  first  necessary  to  consider  the  case  of  diso'ete  detection  of  known  signals  in  additive 
white  Gaussian  noise  (AWGN)  [Ref.  12].  The  power  spectrum  density  of  white  Gaussian 
noise  is  a  constant  over  the  complete  frequency  range.  Since  we  consider  both  negative  and 
positive  frequencies  its  magnitude  is  Nojl.  As  the  probability  density  function  of  a  discrete 
random  process  is  related  to  its  autocorrelation  by  a  Z-transform,  the  corresponding 
autocorrelation  function  for  this  probability  density  function  is  a  delta  function  given  by 
N„/2  •  6(/i)  [Ref.  12]. 

It  is  presumed  that  the  aim  of  this  exercise  is  to  detect  if  a  signal  is  present  or  not.  The 
two  hypotheses  are  given  by: 

Ho:  r[h\  =  n[n]  Noise  only. 

Hi:  r[n]  =  •Si[n]  +  n[n]  Signal  S|  plus  noise. 

In  terms  of  matrix  formulation,  these  are  written  as 

Hq:  r  =  «  Noise  only, 

H,:  r  =  Sj  +  n  Signal  Sj  plus  noise, 

where  the  bold  lower-case  letters  represent  vectors.  Bold  upper-case  letters  are  used  to 
indicate  that  the  symbol  in  question  represents  a  matrix. 
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This  is  a  binary  detection  scenario  as  shown  in  Figure  3.3. 


Ho 

^  Quantity  Threshold 

H. 

t 

Noise[n] 

Figure  3.3:  Binary  detection 

As  shown  in  the  figure,  the  optimum  receiver  computes  the  detection  quantity,  usually 
the  likelihood  ratio,  and  compares  it  to  a  threshold  to  decide  which  of  the  hypotheses 
provides  the  best  answer.  As  we  are  dealing  with  discrete  variables,  the  real  life  input  (i.e., 
analog  input)  is  sampled.  As  discussed  in  Chapter  n,  the  sampling  process  changes  the 
characteristics  of  the  white  noise  to  one  that  is  bandlimited.  As  a  result  of  this  process,  the 
correlation  matrix  of  the  noise  is  no  longer  a  delta  but  rather  takes  the  form  of  a  digital  sine 
function. 

In  order  to  perform  the  detection  using  the  standard  likelihood  ratio  detector  we 
require  that  the  noise  samples  be  uncorrelated.  When  the  sampling  frequency  can  be 
chosen,  one  way  of  ensuring  that  the  noise  samples  are  uncorrelated  is  to  use  a  sampling 
procedme  where  the  sampling  interval  corresponds  to  the  distance  from  the  origin  to  the 
first  zero  crossing  of  the  noise  correlation  function.  Since  the  noise  samples  are 
uncorrelated  and  Gaussian,  they  are  independent.  When  data  is  provided  with  a  set 
sampling  frequency,  the  sampling  interval  can  still  be  set  by  interpolation  and  resampling. 
The  Gaussian  nature  of  the  samples  allows  completely  statistical  description  using  the  first 
two  moments. 
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Therefore,  if  the  sampling  interval  is  chosen  to  occur  at  the  first  zero  crossing  of  the 
correlation  matrix  all  noise  samples  are  uncorrelated  and  the  joint  probability  function 
between  the  signal  and  the  noise  component  of  the  received  discrete  data  can  be  expressed 
as  a  product  of  two  independent  terms.  This  then  allows  the  use  of  a  standard  detector. 
Detectors  designed  to  perform  under  the  assumption  of  white  noise  are  comprehensively 
addressed  in  [Refs.  10, 1 1, 12]  and  are  not  repeated  here. 

To  devise  the  likelihood  ratio  detector,  the  mean  and  variance  under  each  hypothesis 
must  be  calculated.  Since  the  signals  are  known,  the  variance  under  each  hypothesis  is  due 
only  to  the  noise  and  is  therefore  the  same.  This  is  also  the  basis  of  the  whitening  property 
of  the  prediction  error  filter  addressed  later  on  in  this  ch^ter. 

B.  DETECTION  IN  COLORED  NOISE 

When  the  power  spectral  density  of  the  noise  is  not  a  constant  across  the  frequency 
range  it  is  called  colored  noise.  This  means  that  the  correlation  matrix  of  the  noise  is  not  a 
delta  function  and  the  noise  samples  are  correlated.  In  many  signal  processing  applications, 
it  is  assumed  that  the  noise  samples  are  imcorrelated.  Therefore  it  is  often  necessary  to 
transform  a  vector  of  observations  with  correlated  noise  samples  to  one  in  which  they  are 
uncorrelated.  As  a  result  of  such  a  transformation,  the  new  correlation  matrix  of  the 
observations  is  diagonal. 

As  is  the  case  for  white  noise,  the  detection  of  a  signal  corrupted  by  colored  noise  is 
equivalent  to  determining  which  hypothesis  is  true.  The  hypotheses  are 

Hq:  r[h\  =  n^[h\  Colored  noise  only. 

Hi!  r[n]  =  Ji[/i]  +  np[n]  Signals,  plus  colored  noise, 

where  r  is  the  received  signal,  ^2  is  the  known  signal  and  is  the  Gaussian  colored 

noise  of  zero  mean,  having  a  correlation  matrix  denoted  by  .  Baye’s  detection  formula¬ 
tion  provides  the  general  equation  (i.e..  Likelihood  ratio  function  and  threshold),  which 
when  using  different  detail  information  leads  to  a  host  of  different  detection  criteria  (i.e.. 


26 


maximum  likelihood,  min-max,  maximum  a  posteriori,  etc.).  Basic  detection  theory  for 
Gaussian  colored  noise  is  addressed  in  the  following  paragraph  using  the  maximrim  a  pos¬ 
teriori  criterion. 


1.  Maximum  A  Posteriori  Probability  (MAP)  Criterion 

The  detection  technique  whereas  a  decision  must  be  made  as  to  which  one  of  binary 
or  m-ary  hypotheses  is  to  be  selected,  is  called  hypothesis  testing.  In  the  case  of  target 
detection  the  two  hypotheses  are  generally  labelled  as: 

Hq:  Target  not  present. 

Hi!  Target  present. 

For  the  detection  process,  a  decision  has  to  be  made  on  the  basis  of  observed  random 
variables  consisting  of  one  or  more  observations.  Various  rules  can  be  used  for  the 
detection  and  are  derived  by  maximizing  a  measure  of  performance.  The  maximum  a 
posteriori  (MAP)  decision  rule  usually  leads  to  a  partition  of  the  decision  space  into  two 
regions  Ro  and  Rj.  To  illustrate  the  MAP  decision  rule,  we  assume  that  a  choice  must  be 
made  between  Ho  and  Hi  based  on  one  observation  of  the  random  variable  Y  where  y  is  its 
actual  value  at  the  time  of  observation.  The  probability  density  functions  of  Y  related  to 
each  hypothesis  are  known  and  are  /y|//i(3'l^l)  •  Th®  ^  posteriori 

probability  that  Hj  is  the  correct  hypothesis  given  the  observation  y  is  denoted  by  P(//,|y) , 

for  I  =  0, 1 .  These  are  conditional  probability  density  functions  that  are  conditioned  on 
the  observation  [Ref.  13].  The  a  priori  probabilities  will  be  denoted  as  P(Hf),  for 
i  =  0, 1 .  These  factors  denote  the  known  probabilities  of  each  symbol  being  transmitted. 

With  Bayes  rule,  P(/f,|y)  can  be  written  as 


friy) 


(3.1) 


It  is  possible  to  maximize  the  probability  of  correct  detection  by  deciding  in  favor  of 
the  hypothesis  with  the  largest  a  posteriori  probability; 
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Choose  Ho  if  P(//o|}’)  >  j}))  or 

choose  Hi  if  /’(//jly)  >  P(7^q|>>). 
These  decision  rules  can  also  be  written  as 

"o 

With  Equation  (3.1)  this  becomes: 

Hi 

/K|//.(y|^i)  >  pjH^) 
/y|//o(y|^o)  <  ^(^i) 

Ho 


where  the  left-hand  expression  is  called  the  likelihood  ratio,  and  the  right-hand  term  the 
decision  threshold.  The  MAP  decision  rule  then  consists  of  comparing  the  likelihood  ratio 
to  the  detection  threshold.  Since  this  decision  is  based  on  the  conditional  probabilities, 
four  possibilities  exists  in  making  the  determination: 

^  when  Ho  is  true:  No  target,  no  detection. 

:  Choose  Hi  when  Hi  is  true:  Target  detected. 

•  Choose  Ho  when  Hi  is  true:  Target  not  detected. 

P{D^  |//q)  :  Choose  Hi  when  Ho  is  true:  False  alarm, 
where  the  P(Z),  |//p  represents  the  conditional  probability  that  decision  D,-  was  talfP'n 
given  that  is  the  correct  hypothesis. 


The  first  two  probabilities  denote  correct  decisions  with  identified  as  the 

probability  of  detection  (Pp).  P(Di|//q)  is  called  iht  false  alarm  probability  (Pfa)  and 
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P(Di^Hq)  probability  of  a  miss  (Pm).  Each  of  these  regions  is  shown  in  Figure  3.4,  where 
/l(y)  denotes I //j)  and  /qCj)  denotes /y||^^(yj//o) . 


Figure  3.4:  Detection  regions:  Gaussian  probability  density  functions 


In  the  binary  detection  problem,  the  receiver  makes  its  demsion  by  observing  a  number 
of  m  samples  and  then  arrives  at  a  decision  by  comparing  the  information  extracted  to  a 
decision  threshold.  The  Gaussian  colored  noise  case  is  different  from  the  Gaussian  white 
noise  problem  in  that  the  samples  of  the  colored  noise  are  correlated.  The  noise  has  a  non¬ 
constant  power  spectral  density  function  S„(ci))  .In  thatinterval  it  isnot  constant  and  varies 
according  to  the  ambient  noise  affecting  the  observations.  The  MAP  decision  rule  is  then: 

fT\H,(yWi)  >  p(Hq) 

«0 
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Since  is  an  m  x  m  Toeplitz  matrix,  it  can  be  decomposed  into  m  distinct  eigenvalues 

^2, ...» A„  and  m  orthonormal  eigenvectors  e|,  «2> ...» ««  with  the  following  proper¬ 
ties: 
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1,  *  =  / 


(3.11) 


«/  = 


0,  k*l 


Rx^i  ~ 


If  we  group  the  eigenvectors  in  a  matrix  as 

I  I  ...  l' 
•••  ’ 
I  I  ...  I 


t 

E  = 


and  the  eigenvalues  as. 


Aij  0  0  ...  0 

0  A.2  0  ...  0 

0  0  ...  0 


(3.12) 

(3.13) 


(3.14) 


(3.15) 


|_0  0  0  ... 

then 

-  I  (3.16) 

and 

E*^R^E  =  D^.  (3.17) 

Substituting  Equations  (3.11),  (3.12),  (3.13),  (3.16),  and,  (3.17)  into  Equation  (3.10) 
results  in 
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where  the  apostrophe  indicates  that  an  orthogonalization  transformation  is  required  given 

by; 

=  (3.19) 

and 

yo.k  =  ^k^y-  (3.20) 

These  transformations  change  the  colored  noise  detection  problem  to  an  equivalent 
white  noise  detection  problem.  The  decision  rule  shown  in  Equation  (3.18),  is  similar  to  the 
white  noise  scheme  with  one  exception.  In  the  case  of  white  noise,  the  noise  components 
of  the  observations  where  assumed  to  have  equal  variance.  In  the  case  of  colored  noise,  the 
orthogonal  components  of  the  noise  do  not  have  equal  variance  and  this  difference  requires 
some  normalization  which  is  accomplished  in  the  decision  rule  by  normalizing  with  the 
eigenvalue  associated  with  each  eigenvector. 

2.  Basis  of  Whitening  Transformations 

In  the  following  sections,  several  whitening  transformations  are  presented.  It  is  shown 
that  by  using  various  characteristics  of  a  generic  process  jc,  i.e.,  power  spectral  density  or 
correlation  matrix,  that  we  can  whiten  the  power  spectral  density  of  x.  This  transformation 
is  identified  as 

y  -  Ax,  (3.21) 

where  A  is  the  whitening  linear  transformation.  If  x  is  colored  noise  only  (i.e.,  Hq  is  true), 
the  transformation  can  be  written  as 

»c  •  (3.22) 
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where  and  are  colored  and  white  noise  respectively.  If  x  consists  of  signal  plus 
noise  we  get 

y  =  A„(s  +  ii^),  (3.23) 

or 

y  =  A„s  +  «^.  (3.24) 

The  sequence  y  now  consists  of  a  transformed  version  of  the  signal  and  white  noise. 

Therefore  the  goal  becomes  finding  a  matrix  such  that  the  noise  is  whitened. 

There  are  two  fundamental  ways  in  which  a  random  vector  with  correlated  components  can 
be  transformed  into  one  with  imcorrelated  components.  Because  the  correlation  matrix  of 
such  a  process  is  necessarily  diagonal,  this  operation  is  often  called  diagonalization  of  the 
correlation  matrix  [Ref.  14]. 

The  first  method  is  through  the  use  of  an  eigenvector  factorization  of  the  correlation 
matrix.  The  eigenvector  factorization  leads  to  a  diagonal  matrix  of  eigenvalues  a  matrix 
of  eigenvectors.  We  note  that  by  definition  the  eigenvectors  are  orthogonal.  The  second 
method  uses  decompositions  of  the  correlation,  covariance,  or  data  matrix  by  triangular 
factorization. 

Since  we  are  dealing  with  optimal  discrete  time  detection,  the  next  section  starts  with 
the  derivation  of  the  discrete  matched  filter.  The  general  case  where  the  correlation  matrix 
of  the  noise  and  the  signal  are  known  is  derived.  It  is  shown  that  a  solution  to  the  matched 
filter  problem  is  an  eigenvalu^eigenvector  problem. 

C.  DISCRETE  TIME  MATCHED  FILTER 

1.  Deterministic  Signal  in  Colored  Noise 

Let  x[n\  be  a  finite  length  disCTete  vector  of  data  to  be  processed.  It  consists  of  a 
deterministic  signal  s[n]  corrupted  by  additive  colored  noise  n^[h\ : 

x[n]  =  s[n]+n^[n].  (3.25) 
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The  goal  is  to  design  a  Finite  Impulse  Response  (FIR)  filter  to  optimally  filter  (in  the 
sens  of  the  matched  filter)  the  data.  We  want  the  filter  to  pass  the  signal  while  rejecting  the 
noise.  The  block  diagram  is  shown  in  Figure  3.5, 


jf[n] 


y[n] 


Figure  3.5:  Matched  filter 


where  h  [n]  represents  the  filter  impulse  response  and  y  [/i]  is  the  result  of  the  filtering. 
The  output  is 

y[n]  =  +  )’„["].  (3.26) 

where  the  symbol  0  identifies  a  convolution,  y^[/i]  is  the  response  due  to  the  signal  and 

y„[n\  is  the  response  due  to  the  noise.  Obviously  to  be  successful  we  must  maximize  the 

former  while  trying  to  minimize  the  latter.  Therefore,  the  filter  operation  is  to  maximize 
the  ou^ut  signal  to  noise  ratio. 

An  optimal  solution  can  be  obtained  by  the  following  derivation  maximizing  the  SNR 
[Ref.  14].  By  definition. 


SNR  =  -  . 


(3.27) 


Since  the  duration  of  the  signal  is  of  length  /* ,  it  is  assumed  that  the  filter  is  of  a  length 
sufficient  to  allows  it  to  process  the  complete  signal.  We  define  the  following  four  vectors: 
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(3.28) 


(3.29) 


(3.30) 


h[no] 

h[ni] 


h[np_i] 


From  Equation  (3.26)  and  the  definition  of  convolution, 


(3.31) 


}’[«/> -iJ  =  -A:]  =  h  X  , 

t  =  0 


(3.32) 


3’i[”/>-i]  =  Y,  h[k]s[np_i-k]  =  h  5 , 
it  =  o 


(3.33) 


yni^p-i]  =  Y,h[k]n^[np_i-k]  =  h  h^, 
Jt  =  o 


(3.34) 


where  the  ~  indicates  time  reversal  of  the  vector.  With  the  last  two  expressions,  the  SNR 
given  by  Equation  (3.27)  can  be  expressed  as 


£[(A^«,)  («,  li)] 


(3.35) 


h*  s*s  h 

T 

h*  R,h 


(3.36) 
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In  order  to  perform  the  optimal  detection,  we  must  find  the  impulse  response  h  which 

T 

maximizes  the  SNR.  If  we  set  h*  R^Ji  =  1 ,  we  can  concentrate  on  maximiring  the 
numerator.  This  does  not  affect  our  answer  since  this  condition  only  serves  to  normalize 
the  value  of  the  impulse  response.  Therefore,  the  signal-to-noise  ratio  becomes 

SNR  =  h*^rfh.  (3.37) 

The  vector  s  is  the  signal,  so  our  task  is  to  find  the  correct  h .  We  can  do  this  by  the 
use  of  the  Lagrange  multiplier  [Ref.  14].  We  want  to  maximize  the  quantity 

L  =  h*^s*fh  +  Ml  -  h*^R^h) .  (3.38) 

In  this  case,  since  the  constraint  is  real,  X  is  a  real  Lagrange  multiplier.  To  find  the  mini¬ 
mum  we  must  then  find  the  vector  h  for  which  the  gradient  of  L  with  respect  to  is 
zero; 

=  s*f  h-kR„h  =  0.  (3.39) 

By  restating  the  equality  we  get 

is*f)h  =  XR„h,  (3.40) 

which  is  a  generalized  eigenvalue  problem.  It  follows  from  this  that  h  must  be  a  general¬ 
ized  eigenvector  of  Equation  (3.40).  Since  ft  is  a  solution  to  Equation  (3.40),  and  by  using 

h*^  Rjfi  =  1 ,  then 

h*^s*fh  =  Xh*^R^h  =  X.  (3.41) 

This  last  result  shows  that  the  SNR  related  to  the  generalized  eigenvector  ft  is  equal  to  the 
eigenvalue  X  corresponding  to  that  eigenvector. 

It  can  be  shown  that  the  matrix  s*s  is  of  unit  rank  [Ref.  14],  therefore  it  is  possible 
to  find  P-1  linearly  independent  vectors  each  of  which  is  orthogonal  to  the  vector  s*  of 
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length  P.  Since  each  of  these  independent  vectors  are  also  orthogonal  to  s* ,  each  must  be 
an  eigenvector  with  a  corresponding  eigenvalue  of  zero.  Therefore,  the  resulting  SNR 
related  to  these  eigenvectors  is  also  zero.  We  will  show  next  that  the  eigenvector  related  to 

the  largest  eigenvalue  is  proportional  to  R~^s*[Ref.  14;  p.  243].  This  is  done  by 

transforming  the  generalized  eigen-decomposition  problem  of  Equation  (3.40)  to  an 
ordinary  eigen-decomposition  problem. 

the  matrix  is  expanded  in  terms  of  its  eigenvalue  and  eigenvector  matrices 

(3.42) 

and  we  let 

*  =  (3.43) 

and 

(3.44) 

Inserting  the  last  three  Equations  into  Equation  (3.40)  results  in 

(«*^)g  =  Xg  (3.45) 

which  is  an  ordinary  eigenvalue  problem.  It  can  also  be  demonstrated  that  the  matrix 

is  of  umt  rank  and  that  it  also  has  only  one  eigenvector  with  a  corresponding  eigen¬ 
value  not  equal  to  zero.  This  eigenvector  must  be  proportional  to  t .  Therefore 


g°^t. 

(3.46) 

Inserting  Equations  (3.43)  and  (3.44)  into  (3.46)  yields 

(3.47) 

which  when  solved  for  h  results  in 

(3.48) 

—1  T 

Since  E  is  orthonormal,  E  =  E*  ,  and 
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I,<xE,D-Je/s\  (3.49) 

or 

hocR-^^s*  .  (3.50) 

-1/^* 

We  can  also  show  that  s  is  proportional  to  the  correct  eigenvector  by  inserting  it 

in  Equation  (3.40)  which  then  results  in 

=  u* .  (3,5i) 


T  “1 

Since  the  equality  stands,  this  shows  that  our  assumption  is  correct  and  that  s  s*  is  a 
scalar  equal  to  the  eigenvalue  X  as  shown  in  the  next  equation: 

s*ifR„^5*)  =  s*K.  (3.52) 

***7  —I  T  1 

Because  is  Toeplitz,  s  R'^  s*  can  be  rewritten  as  s*  If”  s .  Therefore,  we  can  write 
the  maximum  SNR  as 

~  ^MAX  -  (3.53) 

We  now  know  that 

h  =  aR'^^s* ,  (3.54) 

where  a  is  the  normalizing  factor  that  remains  to  be  found.  Inserting  Equation  (3.54)  in 
the  condition  h*  R^Ji  =  1  and  solving  for  a  gives 


Using  Equations  (3.54)  and  (3.55)  the  optimal  filter  is  then  defined  as 


(3.55) 


(3.56) 
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This  solution  can  be  interpreted  as  applying  a  whitening  transformation  to  the 
sequence  and  using  the  reversed  replica  of  the  transformed  known  signal  to  achieve  optimal 
detection. 

If  the  noise  is  white,  its  correlation  matrix  is  a  diagonal  matrix  with  diagonal  elements 


.  In  this  case  Equation  (3.40)  becomes 

(s*f)h  =  Xo^h  ,  (3.57) 

which  is  an  ordinary  eigen-decomposition  problem. 

2.  Random  Signal  in  Colored  Noise 

The  derivation  of  the  previous  section  for  a  deterministic  signal  can  be  extended  to  that 
of  random  signals  [Ref.  14].  Let  the  signal  5[/i]  and  colored  noise  n^[n]  be  random 
processes.  As  before,  it  is  assumed  that  the  signal  and  the  noise  are  independent,  that  we 
have  P  samples,  and  that  the  signal  is  present  in  all  samples.  Starting  this  time  from  the 
signal  to  noise  ratio  for  two  random  processes  [Ref.  14]: 


SNR  = 


(3.58) 


where 


(3.59) 


(3.60) 


and,  and  are  the  correlations  matrices  of  the  signal  and  the  noise  respectively.  As 

before,  the  magnitude  of  h  is  constrained  so  that  =  1 .  The  SNR  relation 

becomes 


SNR  =  £{|y,[np]|^}  =  h*^R^h. 


(3.61) 
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Using  the  L-agrange  multiplier  again,  to  maximize  the  quantity  L ; 

Z,  =  »•'■«,»  + Ml (3.62) 

To  find  the  minimum  we  must  then  find  the  vector  h  for  which  the  gradient  of  L  with 
respect  to  A  is  zero; 

=  R^h  -  XR„h  =  0.  (3.63) 

By  restating  the  equality  we  get 

R,h  =  \R„h ,  (3.64) 

which  is  a  generalized  eigenvalue  problem  which  h  must  satisfy.  From  this  it  can  be  con¬ 
cluded  that  as  is  seen  in  the  deterministic  case,  the  matched  filter  corresponds  to  the  gener¬ 
alized  eigenvector  (h )  with  the  largest  eigenvalue.  Since  R^  cannot  be  decomposed  into 
its  components  when  j[/j]  is  random,  a  simple  expression  for  h  cannot  be  given  analyti¬ 
cally  and  the  generalized  eigenvalue  problem  must  be  solved  numerically  to  find  the  gen¬ 
eralized  eigenvector  relating  to  the  largest  eigenvalue. 

D.  INVERSE  FILTER 

In  detection  problems,  the  goal  is  to  detect  the  presence  or  absence  of  the  signal.  In 
underwater  acoustic  j^plications,  the  data  consists  of  the  signal(s)  of  interest  which  must 
be  detected,  and  the  sum  total  of  all  transformations  and  noises  added  by  the  ocean  and 
the  transducer.  These  effects  can  be  viewed  as  a  channel  and  can  be  modeled  as  a  linear 
system.  Maximization  of  the  probability  of  detection  of  the  signal  requires  an  inverse 
filtering  operation.  If  a  perfect  inversion  process  can  be  designed,  the  resultant  signal  from 
the  inverse  system  will  recover  the  original  signal.  In  underwater  acoustics,  the  major 
effects  are  due  to  the  colored  noise  and  the  sea  water  attenuation.  Therefore,  the  solution  to 
the  detection  problem  is  to  design  a  corrective  system  which  when  applied  to  the  data, 
minimizes  the  effects  of  the  colored  noise  and  yield  a  replica  of  the  original  signal 
embedded  in  white  (residual)  noise.  Then  a  standard  detector  can  be  used  to  obtain  optimal 
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detection.  In  linear  systems  theory  [Ref.  15],  this  corrective  system  is  called  an  inverse 
system.  In  communication  systems,  it  is  called  an  equalizer.  This  inverse  system  will  have 
a  frequency  response  which  is  the  reciprocal  of  the  distortions  resulting  from  the  channel. 
This  results  in  the  use  of  Least-squares  filter  design  methods.  In  our  case,  the  objective  is 
to  remove  the  effects  of  attenuation  and  of  colored  noise  sources  to  ensure  that  the  power 
spectral  density  of  the  additive  noise  corrupting  the  desired  signal  is  white. 

A  colored  noise  process  n^[n\  can  be  represented  by  a  white  noise  process  ny\n\ 
transformed  by  a  linear  filter  that  is  causally  invertible  [Ref.  16]  as  shown  in  Figure  3.6. 


Figure  3.6:  Colored  noise  created  by  linear  invertible  filter  G[e'^^] . 


(3.65) 


(3.66) 


SM  = 


(3.67) 


If  the  noise  is  to  be  whitened  and  the  correct  power  spectral  density  (a^  )  is  achieved. 


the  inverse  filtering  operation  must  take  place  as  shown  in  Figure  3.7  where  the  noise 
n^[n\  is  whitened  by  the  transformation 
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Figure  3.7:  Colored  noise  whitened  by  inverse  filter  . 


If  it  is  assumed  that  the  power  spectral  density  of  the  colored  noise  S’^-Cto)  is  known, 

the  whitening  filter  H  can  be  designed.  If  a  signal  s[n]  is  present  in  the  colored  noise, 
the  transformation  wUl  also  affect  it.  However,  because  it  is  a  linear  transformation,  the 
frequency  of  the  signal  is  not  affected.  Since  only  the  power  spectral  density  of  the  colored 
noise  is  used  in  the  design  of  the  filter,  the  ratio  of  the  spectral  heights  of  the  signal  and 
noise  at  the  signal  frequency  is  preserved  and  the  transformed  noise  is  whitened. 

E.  WHITENING  BY  UNITARY  TRANSFORMATION 

1.  Eigenvector  Factorization 

Since  the  eigenvectors  of  the  correlation  matrix  for  the  random  vector  x  are 
orthogonal,  they  can  be  used  to  perform  a  unitary  transformation  which  whitens,  hence 
orthogonalizes,  the  noise  components  of  x.  To  elaborate  on  this  method,  let  be  the 

correlation  matrix  of  generic  random  vector  x .  R  j.  is  then  related  to  its  eigenvalues  A,/  and 
eigenvectors  in  the  usual  way 

~  *  (3.68) 

which  means  that  the  correlation  matrix  multiplied  by  one  of  its  eigenvector  is  equal  to  the 
same  eigenvector  multiplied  by  the  corresponding  eigenvalue.  Since  (N  X  N)  is  hermi- 

tian  symmetric,  we  can  find  N  orthonormal  eigenvectors  ,  for  /  =  I...N ,  each  of  which 
has  a  corresponding  real  valued  eigenvalue  X/,  for  /  =  I...N  [Ref.  14].  Since  the  eigen- 
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vectors  are  orthonormal,  they  have  the  property 


H  «/  = 


1,  k  =  l 

0,  k*  I 


(3.69) 


Premultiplying  Equation  (3.68)  by  results  in 


(3.70) 


and,  since  A.^  is  a  constant. 


(3.71) 


We  known  that  and  are  orthonormal  from  Equation  (3.69),  then  (3.71)  becomes 


^k*  -  h 


1,  k  =  l 

0,  k*^l 


(3.72) 


If  we  group  the  eigenvectors  column-wise  in  a  matrix. 


I  I  ...  I 

E  —  ^2^2  ...  iff  y 

I  I  ...  I 


then,  since  all  the  eigenvectors  are  orthonormal. 


E*  E  =  /. 


We  can  now  write  Equation  (3.72)  as. 


(3.73) 


(3.74) 


E*  R,E  = 


(3.75) 


where 
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(3.76) 


0  0  ...  0 

0  A*2  0  ...  0 

^k=  0  0X3...0- 

I  I  I  I  I 

0  0  0  ... 

Let  y  =  Ax  be  a  linear  transformation  of  the  vector  x .  Then, 

E[y]  =  E[Ax]  =  AE[x].  (3.77) 

Therefore  the  mean  of  the  new  random  vector  is  a  scaled  version  of  the  mean  of  x . 

For  the  correlation  matrix, 

Ry  =  Eiyy*"^]  =  E[Axx*'^A*'^] ,  (3.78) 

Ry  =  A£[xx*^]A*^,  (3.79) 

and 

Ry  =  AR^*^.  (3.80) 

Correspondingly  the  covariance  matrix  becomes 

Cy  =  AC^*^  .  (3.81) 

If  we  define  the  random  vector  y  as 

y  -  E*^x,  (3.82) 

and  we  know  from  (3.75)  that, 

E*^R^  =  D^,  (3.83) 

then  we  can  say  using  (3.82)  that, 

Ry  =  E*^R^  =  (3.84) 

T 

where  the  matrix  of  eigenvectors  E*  of  1?^  produced  the  linear  transformation  of  x . 
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This  shows  that  the  correlation  matrix  of  the  new  vector  y  is  the  matrix  of  eigenvalues  of 
X  which  by  definition  is  diagonal.  This  is  an  important  transformation  since  the  samples 
of  y  are  uncorrelated.  Since  the  elements  of  are  always  non-negative  and  usually  pos¬ 
itive,  is  at  least  positive  semidefinite  and  usually  positive  definite.  This  transformation 
whitens  the  samples  of  the  random  vector  x . 

If  we  premultiply  Equation  (3.84)  by  E  and  postmultiply  by  E*  ,  given  that 

T  T 

EE*  =  E*  E  =  I  because  of  their  orthonormality,  then 

=  EDyE*^ .  (3.85) 

This  equation  gives  us  the  correlation  matrix  of  x  in  terms  of  its  eigenvectors  and  eigen- 
values.  Since  E  is  orthonormal,  E  =  E*  ,  and 

=  ED^-^E*^.  (3.86) 

Since  is  a  real  element  diagonal  matrix  with  elements  X,- ,  for  i  =  1...N,  is 

also  a  diagonal  matrix  but  with  elements  l/X^ ,  for  i  =  1...N .  Since  the  X,-  of  are 

all  greater  than  zero  when  R^^  is  positive  definite,  the  diagonal  elements  of  also  are 

greater  than  zero.  Therefore,  this  equation  provides  us  with  a  simple  way  to  invert  the  cw- 
relation  matrix. 

2.  Singular  Value  Decomposition 

Egenvector  factorization  works  well  when  the  correct  correlation  matrix  of  the 
process  described  by  x  is  known.  In  practice  however,  the  correlation  matrix  Rj.  must  be 
estimated  from  a  finite  length  of  data  or  by  an  inverse  Fourier  transform  of  the  estimated 
power  spectral  density.  This  may  result  in  a  matrix  which  is  poorly  conditioned  or  singular. 
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Since  is  generally  estimated,  the  eigenvalues  and  eigenvectors  computed  for  the 
matrix  may  lack  precision.  A  better  way  to  compute  the  eigenvalues  and  eigenvectors  of 
Rj.  is  by  using  the  singular  value  decomposition  of  the  data  matrix  X  (KxM)  defined  as 


^Af+l 


■^M  +  2 


(3.87) 


^(^-l)M  +  l  ^iK-l)M  +  2  ••• 


The  singular  value  decomposition  theorem  states  that  any  A"  x  M  matrix  X  can  be  decom¬ 
posed  into  the  following  product  of  matrices 


X  =  ULV*  , 


(3.88) 


where  U  is  KxK  orthonormal  matrix  containing  left  singular  vectors  arranged  column¬ 


wise, 


I  I  ...  I 

U  =  «j  112  ...  11^  , 

I  I  ...  I 


(3.89) 


V  is  an  M  X  M  orthonormal  matrix  of  right  singular  vectors. 


^  =  Vj  V2  ... 

I  I  ...  I 


(3.90) 


and  L  is  a  AT  X  Af  matrix  of  non-negative  real  singular  values 
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Oi  0 

0  ^2 


0  0 
0  0 


0  0 

This  last  matrix  is  written  for  K>M.  The  a,-  are  in  decreasing  order  and  may  equal  zero 
for  the  larger  values  of  i .  Generally  there  are  r  nonzero  singular  values  o  where  r  is  the 
rank  of  the  matrix  X  (number  of  independent  columns).  If  /C  =  Af  we  have  a  Hiagonai 
matrix.  When  K  <M  the  matrix  L  has  coliunns  of  zeros  rather  than  having  rows  of  zeros, 

and  the  rank  r  of  AT  is  equal  to  its  number  of  independent  rows.  Whether  the  Hatp  is  real 
or  complex,  the  singular  values  are  always  real  and  greater  than  or  equal  to  zero. 

The  equivalence  of  this  technique  to  the  eigenvalue-eigenvector  factorization  method 
is  given  as  follows.  By  definition, 

=  (3.92) 

and  using  Equation  (3.88) 


(3.93) 


T 

U  is  orthonormal,  that  is  U*  U 


I ,  allowing  the  last  step  in  Equation  (3.93). 


This  represents  the  factorization  of  in  the  same  form  as  the  standard  eigenvalue/eigen- 
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vector  factorization. 


=  EDyE*^ ,  (3.94) 

or  by  postmultiplying  both  sides  by  iS , 

R^E  =  EDy^ .  (3.95) 

Since  this  decomposition  is  iinique  we  can  make  the  following  statements;  The  matrix  of 
right  singular  vectors  V  is  the  same  as  eigenvector  matrix  E,  and,  the  singular  values  of 

the  matrix  E  are  the  square  roots  of  the  eigenvalues  if  the  number  of  samples  are 
accounted  for. 

Therefore  the  singular  value  decomposition  can  be  used  as  an  alternative  way  of 
finding  the  eigenvalues  and  eigenvectors  of  the  correlation  matrix  of  a  function  without 
having  to  estimate  its  correlation  matrix.  Because  it  provides  a  better  evaluation  of  the 
eigenvalues  and  eigenvectors,  the  SVD  is  the  method  of  choice  for  solving  most  linear- 
least-square  problems.  It  provides  a  solution  when  other  techniques  such  as  Gaussian 
elimination  or  LU  decomposition  fail  as  a  result  of  singular  «■  close  to  singular  matrices. 
The  SVD  can  be  done  with  a  singular  matrix,  and  provides  a  solution  that  is  “almost” 
unique.  Of  course  to  be  able  to  perform  the  SVD  we  first  must  have  the  data  available  to 
form  the  matrix,  where  the  data  must  be  signal  free  (i.e.,  noise  only). 

3.  Mahalanobis  Transformation 

The  Mahalanobis  transformation  is  another  transformaticMi  which  can  be  used  for 
whitening.  It  is  also  based  on  the  eigendecomposition  of  the  noise  correlation  matrix  and  is 
defined  as  follows  [Ref.  14:  p.  247]: 

y  =  X  =  .  (3.96) 

This  transformation  results  in  a  diagonal  correlation  matrix  of  the  random  variable  y 
as  shown 
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(3.97) 


Inserting  Equation  (3.83)  results  in 

Ry  =  EDj^'  .  (3.98) 

and 

Ry  =  EIE*^  =  /,  (3.99) 

proving  that  the  transformation  shown  in  Equation  (3.96)  results  in  a  unit  variance  white 
noise  random  process. 

F.  WHITENING  BY  TRIANGULAR  FACTORIZATION 

As  previously  discussed,  transformations  that  whiten  colored  noise  must  result  in  a 
correlation  matrix  which  is  diagonal.  Another  way  of  achieving  this  effect  is  by  triangular 
factorization,  also  known  as  triangularization.  The  term  tiiangularization  refers  to  the  form 
of  the  matrix  of  orthogonal  vectors  which  is  used  to  orthogonalize  the  colored  noise 
component  of  the  received  samples.  Otherwise  this  technique  is  similar  to  whitening  by  a 
unitary  transformation.  This  method  leads  to  a  transformation  that  can  be  interpreted  as 
causal  or  anticausal  linear  filtering  of  the  associated  signal  sequence  depending  on  the  use 
of  the  lower-upper  triangularization  or  the  upper-lower  triangularization. 

1.  Lower-Upper  Factorization  (LDU) 

A  square  matrix  R^^  can  be  expressed  as 

R^  =  LDjU,  (3.100) 

where  L  is  a  unit  lower  triangular  matrix,  meaning  that  the  diagonal  elements  are  equal  to 
1,  and  all  elements  above  the  diagonal  are  equal  to  zero,  Dj^  is  a  diagonal  matrix,  and  V 

is  a  unit  upper  triangular  matrix.  Since  the  correlation  matrix  R^^  is  symmetric,  this  trans¬ 
formation  takes  the  form 
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(3.101) 


=  LDj^L*^. 

The  correlation  of  a  process  can  be  diagonalized  by  the  following  transformation; 

y  =  L’^x.  (3.102) 

To  show  that  this  transformation  diagonalizes  the  correlation  matrix  of  y  we  begin  with 
the  correlation  function  fory, 

Ry  =  E{yy*^}  (3.103) 

using  Equation  (3.102), 

Ry  =  E.  (L^x)iL^x)*^  (3.104) 

=  L^E{xx*^HL^)*^,  (3.105) 

and, 

=  L-^R^iL~^)*^ .  (3.106) 

If  Equation  (3.101)  is  now  inserted  in  Equation  (3.106), 

Ry  =  L-\LDj^L*^)iL~^)*^ .  (3.107) 

Now  L  L  and  L*  (L  )  cancel  out  and  we  are  left  with 

Ry  =  Di^.  (3.108) 

Since  we  know  that  is  a  diagonal  matrix  this  proves  that  the  transformation  does 
indeed  diagonalize  the  correlation  matrix  of  y .  If  is  the  correlation  matrix  of  a  colored 

noise  process,  then  the  transformation  of  Equation  (3.102)  whitens  it 
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When  applied  to  a  sequence  of  consecutive  samples  of  a  random  sequence,  this 
transformation  is  causal  and  results  in  a  series  of  orthogonal  random  variables. 


2.  Upper-Lower  Factorization  (UDL) 

An  alternative  triangular  decomposition  factors  the  correlation  or  covariance  matrix 
into  a  upper-lower  decomposition.  This  decomposition  of  a  correlation  matrix  has  the  form 

=  UDjjU*^,  (3.109) 

where  U  is  unit  upper  triangular  and  the  matrix  Djj  is  diagonal.  This  last  matrix  is  differ¬ 
ent  than  the  matrix  of  the  lower-upper  decomposition. 


If  we  premultiply  the  previous  equation  by  U~^  and  postmultiply  by 


tT-  -1 


iu*‘)  =(l/'^) 


-1 


-1  -1 


(3.110) 


Therefore  keeping  in  mind  in  mind  the  definition  for  ,  it  can  be  shown  that  the  transfor¬ 

mation 

y  =  u~^x  ,  (3.111) 

also  results  in  a  vector  with  orthogonal  components.  Since  U  is  upper  triangular,  this 
transformation  is  anti-causal.  The  matrix  1/  can  be  computed  with  the  LDU  method  in  the 
following  way  [Ref.  17]: 

Taking  the  reversal  of  Equation  (3.109)  results  in 

Rjc  =  UDjjU*^,  (3.112) 

where  the  matrix  U  is  lower  diagonal.  Equation  (3.112)  then  has  the  LDU  form  of  Equa¬ 
tion  (3.1(X)).  Since  this  decomposition  is  unique  it  means  that  the  lower  triangular  matrix 

resulting  LDU  decomposition  of  the  matrix  R^  is  equal  to  U .  The  matrix  U  of  Equation 
(3. 1 1 1)  is  then  found  by  reversing  the  lower  triangular  matrix  achieved  from  Rx . 
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3.  Cholesky  Factorization 

The  Cholesky  decomposition  is  a  special  case  of  the  LDU  decomposition  [Ref.  16].  If 
can  be  expressed  as 

R^^LDjL*^,  (3.113) 

and,  the  diagonal  elements  of  are  positive,  the  matrix  R^  can  then  be  further 
expressed  as 

Rx  ~  ^c^C*  ’  (3.114) 

where  is  lower  triangular.  Therefore  the  correlation  matrix  of  a  process  can  be  diago¬ 
nalized  by  the  transformation 

y  =  Lc^x,  (3.115) 

which  is  causal  since  Lq  is  lower  triangular. 

In  the  same  way,  starting  from  the  UDL  decomposition  we  can  also  generate  a 
Cholesky  factorization  which  results  in  the  factorization  of  into  an  upper  triangular 
matrix  followed  by  a  lower  triangular  matrix, 

Rx  -  UqU(^*  ,  (3.116) 

where 

Uf.  =  (Ic)*’'.  (3.117) 

and  which  leads  to  the  anti-causal  transformation 

y  =  Uc^x.  (3.118) 

Since  both  transformations  at  Equation  (3.115)  and  (3.118)  are  special  cases  of  the 
LDU  factorization  it  can  easily  be  shown  that  they  also  whiten  a  colored  noise  sequence. 
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4.  QR  Factorization 

The  Cholesky  and  LDU  decompositions  provide  a  means  by  which  the  normal 
equation  can  be  solved  by  performing  square  root  decompositions  of  the  autocorrelation 
matrix  .  The  QR  factorization  decomposes  a  data  matrix  X  into  an  orthogonal  matrix  Q 

and  an  upper  triangular  matrix  R .  Since  this  decomposition  uses  the  data  matrix  directly, 
it  provides  better  numerical  accuracy  than  the  square  root  methods  which  are  based  on  the 

decomposition  of  the  autocorrelation  matrix  R^^  which  is  usually  estimated  from  a  finite 
length  data  [Ref.  17]. 

Let’s  consider  the  following  derivation  from  [Ref.  14].  To  begin,  we  consider  the 
KxM  data  matrix 


X  = 


I  I  ...  I 

Xi  X2  ...  Xj^  , 

I  I  ...  I 


(3.119) 


where  x^  =  [jt.[i]  a:, [2]  ...  .^,.[/:]]^, for  i  =  l...M,andthematrixisoffullrank, 

meaning  that  the  numbers  of  rows  A"  is  at  least  equal  to  the  number  of  columns  M  and  that 
these  latter  are  independent.  Since  the  columns  are  independent,  the  Gram-Schmidt 
orthogonalization  procedure  can  be  used  to  derive  a  set  of  M  orthonormal  vectors  >  92  > 

93  .•..»93f  •  For  the  Gram-Schmidt  method,  we  choose  a  normalized  version  of  the  first 
column  as  the  first  orthonormal  vector: 

9'l  =  (3.120) 

and  normalize  it. 


= 


(3.121) 


To  form  the  next  orthonormal  vector,  we  take  the  second  coliunn  and  remove  from  it  the 
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Pll  Pl2  —  Piw| 

0  P22  —  P2Af 

. ,  (3.126) 

0  0  ... 

where  represent  the  transformation  coefficients  that  are  q)plied  to  the  data  vectors  to 
create  the  orthonormal  vectors  qi .  This  equation  can  also  be  written  as 

Q  =  XR^  (3.127) 

If  we  multiply  both  sides  of  Equation  (3.127)  with  R , 
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X  =  QR.  (3.128) 

This  last  equation  now  expresses  the  matrix  AT  as  a  function  of  a  rectangular  matrix  Q 
whose  columns  are  orthonormal,  and  a  matrix  R  which  is  square  upper  triangular. 

To  show  that  the  QR  decomposition  can  be  used  to  obtain  the  triangular 
decomposition  of  an  estimated  correlation  matrix  defined  by 

«x  =  (3.129) 

With  (3.128)  this  becomes 

R^  =  ^R*^Q*^QR  .  (3.130) 

Since  (2  is  an  orthonormal  matrix, 

Q*^Q  =  I.  (3.131) 

Therefore  Equation  (3.130)  becomes 

*x  =  (3.132) 

Rj^  can  also  be  factored  by  LDU  decomposition. 


So  it  follows  that 

R^  =  LDj^L*^ . 

(3.133) 

=  LD^L*\ 

(3.134) 

resulting  in 

and. 

R*^  = 

(3.135) 

R  = 

(3.136) 

If  this  equation  is 

1/2 

solved  for  ,  we  have 
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(3.137) 


where  the  formula  for  Dj^  is  the  previous  one  given.  Then  we  can  do  the  same  transfor¬ 
mation  as  for  LDU, 

y  =  L~^x  .  (3.140) 

Since  the  QR  factorization  directly  uses  the  data  matrix  to  compute  the  triangular 
matrix  ,  it  is  not  subject  to  the  problems  resulting  from  using  an  estimated  correlation 
matrix.  This  makes  this  technique  less  subject  to  errors  resulting  from  poor  conditioning 
and  usually  leads  to  better  numerical  accuracy. 


G.  DIFFERENTIAL  FILTER 

In  Chapter  n,  it  is  shown  that  the  absorption  of  sound  energy  by  sea  water  has  an 
important  effect  on  the  ambient  noise  in  the  ocean.  The  intensity  attenuation  coefficient  is 
defined  as  being  proportional  to  the  square  of  the  frequency  for  most  frequencies.  This 
resulted  in  a  sound  intensity  magnitude  proportional  to  frequency  squared, 
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(3.141) 


and  the  pressure  magnitude. 


(3.142) 


For  the  noise  to  be  white,  it  is  required  that  the  magnitude  of  the  frequency  response 
be  equal  to  a  constant.  The  colored  noise  signal  n^[n\  is  related  to  its  frequency  response 


[®]  by  the  Fourier  transform 


[®1- 


The  frequency  response  is  equal  to 


to 


(3.143) 


(3.144) 


where  y4  is  a  constant.  If  the  random  process  is  differentiated  with  respect  to  time, 
in  the  frequency  domain  this  has  the  effect  of  multiplying  [co]  by  a  factor  of  yto ; 


(3.145) 


__  d  .  271  ycD/ 
dt  to 


(3.146) 


.  271  d  f 

=  ). 

to  dr 


(3.147) 


whereto  =  27t/,  resulting  in 


(3.148) 


This  last  equation  shows  that  the  frequency  spectrum  of  the  time  differentiation  of 
n^[n\  should  have  a  constant  magnitude  equal  to  2nA  and  therefore  has  been  whitened. 
Since  the  differential  operator  is  a  linear  operation,  at  the  frequency  of  a  given  signal,  the 
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ratio  between  signal  power  and  noise  power  should  remain  the  same  and  the  noise  is 
whitened. 

H.  WHITENING  PROPERTY  OF  PREDICTION-ERROR  FILTER 
A  prediction  error  filter  (PEF)  is  shown  in  Figure  3.8. 


Figure  3.8:  Prediction-MTor  filter 


The  PEF  combines  the  features  of  a  predictor  filler  whose  role  it  is  to  estimate  the 
value  of  the  sample  u[n\  based  on  the  M  previous  values  starting  at  a  delay  of  d.  The  error 
between  the  predicted  value  and  the  real  value  is  tirai  calculated  by  the  summer.  Often  a 
variation  of  this  set-up  is  used  which  computes  the  M  weights  of  die  predictor  in  an 
adaptive  fashion.  The  advantage  to  the  adaptive  filter  is  as  its  name  implies  is  that  it  has  the 
ability  to  change  its  behavior  as  the  environment  changes.  For  example  in  the  SONAR 
processing  problem,  the  characteristics  of  the  ambient  nmse  may  change  with  variations  in 
locations,  weather  conditions,  time  of  day,  and,  movmnents  of  noise  generators.  The 
adaptive  filter  will  change  its  weight  to  ensme  that  its  error  is  minimized.  Li  those  cases. 
Figure  3.8  is  modified  as  shown  on  Figure  3.9. 


59 


Figure  3.9:  Adaptive  PEF 


As  shown,  the  errors  are  used  to  update  the  weights  of  the  predictor  to  minimize  the 
difference  between  the  real  signal  and  its  estimate.  This  type  of  filter  relies  on  the 
correlation  characteristics  of  the  components  that  are  being  separated.  In  our  case,  the  two 
components  are  the  noise  and  the  known  signal  that  is  to  be  detected. 

In  the  case  of  white  noise,  noise  samples  with  a  delay  of  one  or  more  between  them 
are  uncorrelated.  Therefore,  the  predictor  can  be  used  to  effectively  estimate  the  values  of 
the  known  signal  if  the  delay  is  chosen  at  a  point  where  the  correlation  function  of  the  signal 
is  at  a  maximum. 

In  the  case  of  colored  noise,  the  noise  and  the  signal  may  have  correlation  function 
which  are  less  distinctive  and  additional  care  must  be  taken  in  the  design  of  the  filter.  The 
correlation  function  of  the  low-pass  Butterworth  noise  for  the  first  25  delays  is  shown  in 
Figure  3.10. 
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Figure  3. 10:  Correlation  function  of  experimental  first  order  low-pass  Butterworth  coltxed 

noise. 

As  can  be  readily  observed,  the  correlation  function  for  colored  noise  is  a  relatively 
smooth  fimction  which  does  not  oscillate  around  zero  at  predetermined  delays.  Diffo'ent 
realizations  of  the  noise  result  in  different  zero  crossovers  but  the  general  form  of  the 
function  is  maintained  and  relatively  strong  correlation  is  present  only  at  short  delays. 
However,  because  the  known  signal  to  be  detected  is  a  periodic  signal,  its  correlation 
function  is  also  periodic  and  as  such  goes  through  zero  at  known  delays.  This  may  allow  a 
delay  to  be  chosen  that  decorrelates  the  signal  while  maximizing  the  correlation  of  tiie 
noise.  If  the  signal  is  decorrelated,  the  predictor  can  provide  a  reasonable  estimate  of  the 
colored  noise.  In  this  case,  the  error  between  the  estimated  noise  value  and  the  real  value 
consists  of  the  signal  and  a  prediction  error.  The  prediction  aror  is  based  on  two  factors; 
the  amount  of  decorrelation  between  the  noise  samples  at  the  chosen  delay,  a 
quantization  error.  The  predictor  relies  on  the  correlation  between  adjacent  samplftig  of  the 
input  process.  Therefore,  when  the  order  of  the  filter  is  increased,  the  correlation  of  the 
adjacent  samples  of  the  input  prcx^ess  is  successively  reduced  until  the  filter  is  of  an  orda 
where  the  output  samples  are  uncorrelated  and  therefore  white  [Ref.  18:  p.  216].  Howeva, 
so  that  the  desired  signal  is  not  whitened,  the  requires  that  the  delay  selected  be  such 
that  the  (x>rrelation  of  the  component  to  be  predicted,  in  our  case  the  colored  noise,  be  as 
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large  as  possible  while  that  of  the  signal  has  a  correlation  as  close  to  zero  as  possible.  Any 
amount  of  correlation  at  that  delay  in  the  desired  signal  causes  the  predictor  weights  to 
converge  to  a  state  where  some  energy  of  the  other  signal  is  present  in  the  estimated 
quantity  and  is  therefore  removed  from  the  “error”  signal.  Therefore,  in  our  problem,  the 
predictor  estimates  the  noise  component  of  the  received  signal.  If  the  desired  signal  is 
uncorrelated  at  the  delay  chosen,  then  the  error  e[n\  consists  in  the  desired  signal  and  a 
white  noise  prediction  error.  The  filter  whitens  the  noise  corrupting  the  signal  if  the  error 
is  chosen  as  the  desired  output.  The  variance  of  the  error  depends  on  the  magnitnHe  of  the 
learning  rate  parameter.  To  a  certain  extent,  the  smaller  the  parameter  is,  the  lower  is  the 
variance  of  the  prediction  error  but  the  filter  weights  of  the  predictor  take  longer  to 
converge  and  are  less  responsive  to  sudden  variations  in  the  acoustic  environment. 
Conversely,  if  the  learning  rate  parameter  is  larger,  the  weights  converge  more  quickly  but 
the  prediction  errors  are  larger  resulting  in  white  noise  of  higher  variance  in  the  predicted 
error  and  therefore  a  lower  output  signal  to  noise  ratio. 
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IV.  COLORED  NOISE  DETECTION  RESULTS 


This  chapter  presents  the  results  achieved  with  the  various  methods  addressed  in 
Chapter  III.  Each  whitening  technique  is  addressed  in  the  order  in  which  it  is  discussed  in 
Chapter  EH  beginning  with  the  matched  filter  and  inverse  filter,  and  closing  with  the 
prewhiteners. 


A.  TEST  PROCEDURES 

A  common  test  sequence  is  used  for  the  testing  of  all  methods  except  the  prediction- 
error  filter  which  uses  a  different  signal.  This  sequence  consists  of  three  sinusoids  of 
different  frequencies  corrupted  by  colored  noise.  The  sinusoids  form  the  known  signal  that 
is  to  be  detected.  They  are  located  in  different  areas  of  the  spectrum  to  allow  a  simultaneous 
view  of  the  effects  at  each  spectral  location.  One  sinusoid  is  at  a  very  low  frequency  where 
the  slope  of  the  noise  is  greater,  one  is  at  higher  frequencies  where  the  noise  spectral  density 
is  almost  flat,  and  one  is  located  at  an  intermediate  frequency.  The  frequencies  of  the 
sinusoids  are  0.05  fs,  0.15  fs,  and,  0.3/s. 

The  noise  is  created  by  passing  a  white  Gaussian  noise  through  a  low-pass  filter.  It  is 
added  to  the  signal  of  interest  and  is  colored  with  a  first  order  low-pass  Butterworth  filter 
with  a  cut-off  frequency  of  0.1  fs.  The  Butterworth  filter  is  an  all-pole  filter.  Its  magnitude- 
square  frequency  response  is  given  by  [Ref.  19] 


l^(co)|^  = 


1 

l  +  (co/©^)^^’ 


(4.1) 


where  N  is  the  filter  order,  ©  the  digital  frequency  and  (s>^  the  3dB  cut-off  digital  fre¬ 
quency.  Therefore,  a  first  order  Butterworth  filter  has  a  frequency  response  proportional  to 

|f/(©)|^oc±,  (4.2) 

© 

which  is  similar  to  the  frequency  response  observed  for  the  ambient  noise  in  the  ocean. 
Since  a  Butterworth  filter  is  used  to  create  the  colored  noise,  the  theoretical  power  spectral 
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density  and  correlation  function  of  the  noise  is  known.  The  power  in  the  signal  is  chosen 
so  that  the  spectral  heights  of  each  of  the  sinusoids  and  the  noise  have  a  ratio  of  approxi¬ 
mately  two  at  the  respective  spectral  locations,  when  the  white  Gaussian  noise  used  as  a 
input  to  the  Butterworth  noise  filter  has  a  variance  of  100.  Therefore,  the  three  sinusoids 
forming  the  signal  have  amplitudes  of  0.475, 0.2,  and,  0.07  respectively.  Different  meth¬ 
ods  require  that  higher  input  signal  to  noise  ratios  be  used.  For  these  results,  the  white 
noise  power,  i.e.,  its  variance,  used  to  create  the  colored  noise,  is  reduced.  This,  in  effect, 
reduces  the  power  of  the  colored  noise  proportionally  across  the  frequency  spectrum. 

Figure  4.1  shows  an  averaged  (100  realizations)  power  spectral  density  of  the  test 
sequence  with  the  variance  of  the  white  Gaussian  noise  set  at  100. 


Figme  4.1:  Normalized  power  spectral  density  of  signal  plus  noise 


The  power  spectral  density  in  the  area  of  Sinusoids  #2  and  #3  are  also  plotted 
separately  on  each  {q>plicable  figure  to  better  illustrate  the  relationship  between  the  signal 
power  and  the  noise  power  at  those  frequencies.  This  is  also  done  in  all  result  plots  where 
relative  power  levels  are  very  different  from  one  part  of  the  frequency  spectrum  to  the 
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other.  All  plots  are  normalized  to  the  largest  peak  since  our  interest  is  in  the  signal  to  noise 
ratios  resulting  from  the  transformations. 

The  power  spectral  density  from  each  processing  scheme  is  shown  using  averaged 
periodograms.  The  periodogram  is  a  well  known  estimate  for  the  power  spectral  density. 
Ten  realizations  of  each  method  are  performed  and  the  ten  results  averaged.  This  is  done 
because  the  pmodogram  is  by  definition  an  unreliable  estimator  due  to  its  large  variance. 
When  the  results  of  ten  periodograms  are  averaged,  the  variance  is  reduced  by  a  factor  of 
ten,  rendering  a  more  consistent  estimate.  In  all  cases  the  estimated  power  spectral  density 
of  the  sequence  before  and  after  whitening  are  shown.  If  the  transformation  effectively 
decorrelates  the  noise  components  of  the  samples,  the  three  sinusoids  afro-  decorrelation 
should  all  have  approximately  the  same  signal  to  noise  ratio  at  their  respective  frequency. 

In  aU  cases  the  processing  consists  of  the  following  steps: 

Step  a:  Add  the  signal  and  colored  noise  to  create  the  test  sequence. 

Step  b:  Derive  the  transformation  vector  or  matrix  depending  on  the  algorithm  being 
tested. 

Step  c:  Operate  on  the  test  sequence  the  applicable  transformation,  and  raeate  a  new 
sequence  in  which  the  noise  has  been  whitened. 

Step  d.  Use  the  fast-Fourier  transform  operation  to  form  a  bank  of  discrete  detectors 
to  attempt  to  locate  the  three  signal  components  in  the  white,  or  whitened  noise.  The  results 
are  computed  as  periodograms. 

Step  e:  Run  the  transformation  of  ten  different  realizations  of  the  e3q)eriment  and 
compute  an  averaged  periodogram. 

In  terms  of  detection,  the  resulting  power  spectral  density,  for  two  signals  in  noise,  Si 
and  $2  in  this  example,  is  interpreted  as  shown  on  Figure  4.2: 
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A  detection  threshold  will  not  be  set  per  se  but  rather  we  will  be  looking  at  how  die 
power  ratio  between  the  signal  and  the  noise  at  the  respective  frequencies  of  the  three 
sinusoids  are  maintained  by  each  transformation. 

1.  Estimates  of  the  Correlation  Matrix 

Many  of  the  methods  that  are  discussed  require  that  the  correlation  function  of  the 
noise  be  known  or  estimated.  We  assume  that  the  power  spectral  density  of  the  ambient 
noise  is  known. 

When  the  power  spectral  density  of  the  noise  is  known,  its  autocorrelation  function  can 
be  estimated  using  the  Wiener-Khintchine  relation  [Ref.  19]. 
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L«t  be  the  power  spectral  density  of  the  wide  sense  stationary  random  vector 

X .  The  autocorrelation  is  related  to  the  power  spectral  density  as  follow: 

oo 

(4.3) 

/=  -OO 

which  is  the  discrete  time  Fourier  transform  of  the  autocorrelation  function.  Since  and 

are  DTFT  pairs,  the  latter  can  be  given  by: 

n 

-n 

It  is  reasonable  to  assume  that  the  noise  power  spectral  density  in  a  given  area  of  ocean 
is  known  with  a  certain  degree  of  accuracy  as  a  result  of  previous  data  and  knowledge  of 
the  present  environmental  conditions.  In  most  of  the  tests,  it  is  assumed  that  the  power 
spectral  density  of  the  noise  is  known.  From  this  power  spectral  density,  the  corresponding 
autocorrelation  function  can  be  derived  using  the  Wiener-Khintchine  relationship. 

2.  Matrix  Decomposition  Methods 

These  methods  are  realized  by  using  block  transformations  to  filter  the  sequence.  This 
means  that  the  resultant  sequence  is  formed  by  one  matrix  multiplication  obtained  via  the 
product  of  a  matrix  of  orthonormal  vectors  and  the  input  sequence.  One  of  the  problems 
with  using  this  technique  is  that  it  is  usually  preferable  to  use  larger  matrices  allowing  more 
input  data.  This  implies  that  the  correlation  matrix  to  be  decomposed  is  also  larger.  The 
processing  resources  required  to  factor  larger  matrices  into  its  required  orthonormal 
constituents  becomes  quickly  prohibitive.  To  attempt  to  mitigate  this  problem,  the 
correlation  matrices  used  in  these  methods  are  limited  to  a  size  of  256  X  256.  This  means 
that  the  multiplications  are  also  performed  in  blocks  of  256  points.  Since  we  want  to  be  able 
to  compare  the  results  achieved  with  all  methods,  and  we  want  representative  results,  all 
power  spectral  densities  shown  are  for  data  lengths  of  1024.  For  matrices  factorization 
methods  limited  by  size,  four  transformed  block  of  data  are  then  joined  to  form  one  1024 
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length  sequence.  It  is  expected  that  the  results  may  be  distorted  since  the  transformed  data 
points  starting  at  n,  /i+256,  n+512,  and,  /i+768,  would  then  all  be  the  results  of  the  same 
orthonormal  vector  multiplying  different  input  sequences.  However,  when  only  four  blocks 
are  used,  this  problem  is  not  observed.  When  smaller  matrices  and  therefore  more  data 
blocks  are  used  in  an  attempt  to  reduce  the  processing  requirements,  a  clear  periodic  pattern 
emerged  based  on  the  number  of  blocks  used.  Therefore,  the  use  of  many  small  block 
transforms  leads  to  unsatisfactory  results  and  is  avoided. 

B.  RESULTS  OF  OPTIMAL  DETECTION 

1.  Discrete  Time  Matched  Filter  for  Colored  Noise 

This  section  addresses  the  results  for  the  discrete  time  matched  filter  for  colored  noise 
as  derived  in  Chapter  m.  The  figures  shown  were  computed  using  the  Matlab  program 
given  in  Appendix  H.  The  results  are  achieved  by  filtering  the  test  sequence  x  with  h  to 
arrive  at  a  sequencey.  The  transfer  function  h  of  the  filter  is  computed  using  Equation  3.56. 


Figure  4.3:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=l(X)). 


Figure  4.3  shows  the  normalized  averaged  power  spectral  density  for  the  input 
sequence  x.  As  stated  previously,  the  three  sinusoids  are  at  fiiequencies  0.05  fs,  0.15  fs,  and, 
0.3  fs.  The  variance  of  the  white  Gaussian  noise  used  to  create  the  colored  noise  is  100.  The 
power  in  each  of  the  three  sinusoids  is  selected  so  that  the  spectral  heights  of  signal  and 
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noise  at  each  frequency  have  a  ratio  of  approximately  two.  Figure  4.3  shows  an  averaged 
periodogram  based  on  ten  realizations  rater  than  100  realizations  as  in  Figure  4.1.  The 
variance  of  the  periodogram  is  evident  when  comparing  Figures  4.1  and  4.3.  In  Figure  4.3, 
Smusoid  #2  and  Sinusoid  #3  are  not  distinguishable  from  the  noise  background.  Sinusoid 
#1  can  be  identified  but  other  peaks  at  lower  frequencies  are  clearly  stronger  which  maifp^ 
its  detection  unlikely  with  the  periodogram 

The  matched  filter  is  usually  used  to  generate  one  number  whose  value  is  compared  to 
a  threshold  to  determine  the  presence  or  absence  of  a  signal.  Because  of  the  equivalence  of 
the  matched  filter  and  the  correlator,  the  output  of  the  matched  filter  can  be  interpreted  as 
a  correlator.  Therefore,  because  we  want  to  compare  aU  results  in  terms  of  averaged 
periodograms,  the  results  of  the  matched  filter  detection  are  shown  in  the  form  of  a 
correlator.  The  averaged  periodogram  is  computed  for  the  sequence  consisting  of  all 
matched  filter  results.  The  next  figure  shows  the  nramalized  averaged  power  spectral 
density  of  the  filtered  sequence  y. 


Figure  4.4:  Normalized  power  spectral  density  of  output  of  matched  filter. 


Figure  4.4  clearly  shows  the  success  of  the  algorithm  in  filtering  the  colored  noise.  The 
strong  relative  power  of  each  sinusoid  to  the  noise  is  well  shown.  All  three  sinusoids  have 
significant  signal  to  noise  ratios,  far  better  than  those  in  Figure  4.3.  The  peaks  are  shaip  and 
precisely  located  at  the  correct  frequencies.  The  effect  of  the  colored  noise  appears  to  have 
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been  almost  completely  eliminated  with  the  exception  of  some  residual  noise  around  each 
sinusoid.  A  detection  threshold  could  be  established  at  a  very  low  level  with  no  chances  of 
false  alarms  or  of  misses.  TTiese  results  show  that  the  matched  filter  fw  colored  noise 
applications  is  very  successful. 

2.  Inverse  Filter 

The  results  achieved  for  the  inverse  filter  derived  in  Chapter  HI  are  shown  here.  These 
results  were  achieved  using  the  Matlab  program  shown  in  Appendix  I.  The  inverse  filter  is 
used  as  a  prewhitener.  The  results  where  achieved  by  running  the  test  sequence  x  through 
the  filter  point  by  point  to  arrive  at  a  sequence  y  consisting  in  mostly  signal.  The  transfer 
function  k  of  the  filter  is  computed  with  Equation  3.66  by  assuming  that  the  power  spectral 
density  of  the  colored  noise  and  therefore  the  transfer  function  reqwnsible  for  the 
colorations  of  the  noise  is  known. 

Figure  4.5  shows  the  normalized  averaged  power  spectral  density  for  the  input 
sequence.  In  this  case,  the  noise  power  is  reduced  by  setting  the  white  Gaussian  noise 
variance  to  25.  Compared  to  the  case  when  the  white  Gaussian  noise  variance  is  set  to  100, 
the  use  of  the  lower  variance  provides  a  gain  of  approximately  6  dB.  This  gain  is  shown 
clearly  in  Figure  4.5  where  the  signal  peaks  are  much  more  clearly  distinguishable  than  in 
Figures  4.1  and  4.3.  In  fact  Sinusoid  #1  could  be  reliably  detected  with  a  simple  detection 
threshold  despite  the  higher  noise  levels  at  lower  frequencies.  Sinusoid  #2  and  #3  however 
have  far  less  power  than  the  low  frequency  noise  and  Sinusoid  #1. 
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Figure  4.5:  power  spectral  density  of  signal  +  colored  noise  (white  Gaussian  noise 

variance=25) 


The  results  of  filtering  the  sequence  x  with  the  inverse  filter  are  shown  at  Figure  4.6: 
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Figure  4.6:  Normalized  power  spectral  density  of  whitened  signal  using  inverse  filter. 


We  observe  that  the  noise  has  been  completely  whitened  as  theoretically  predicted.  It 
can  also  be  seen  that  the  three  components  of  the  signal  can  be  clearly  detected  and  have 
approximately  the  same  relative  SNR  at  their  respective  frequency  locations.  Because  the 
inverse  filter  is  essentially  an  equalizer,  the  lower  power  sinusoids  have  been  “raised”  to 
approximately  the  same  level  as  the  first.  Commensurately,  the  noise  at  low  frequencies  is 
reduced  and  that  at  higher  frequency  is  increased.  The  transformation  does  not  appear  to 
produce  peaks  at  incorrect  frequencies  which  could  lead  to  false  alarms.  Since  the  inverse 
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filter  does  not  require  decomposition  of  large  matrices  or  matrix  multiplications,  the 
processing  requirements  for  this  whitener  are  relatively  minimal. 

3.  Eigenvector  Factorization 

The  eigenvector  factorization  results  are  the  first  of  the  matrix  decomposition  methods 
that  are  shown.  The  results  achieved  for  the  eigenvector  transformation  derived  in  Chapter 
in  are  shown  here.  These  results  were  obtained  using  the  ^4atlab  program  shown  in 
Appendix  J .  The  results  are  achieved  by  separating  the  test  sequence  x  into  four  blocks  and 
multiplying  each  of  these  by  the  orthonormal  matrix  of  eigenvectors  of  the  theoretical  noise 
correlation  matrix  of  size  256  X  256. 

Figure  4.7  shows  the  normalized  power  spectral  density  for  the  input  sequence.  The 
noise  power  is  reduced  by  setting  the  white  Gaussian  noise  variance  to  10.  This  is  done  due 
to  the  poor  performance  of  the  whitener  as  is  shown  on  Figure  4.8.  The  low  noise  power  is 
clearly  apparent  in  Figure  4.7  by  the  clear  imambiguous  signal  peaks  and  the  relative  lack 
of  significant  noise  levels.  In  fact  Sinusoid  #1  can  be  reliably  detected  with  a  simple 
detection  threshold  despite  the  higher  noise  levels  at  lower  frequencies.  We  note  that 
Sinusoid  #2  and  Sinusoid  #3  would  be  also  be  easily  detected  by  common  bandpass  filters 
with  predetermined  detection  thresholds. 


Figure  4.7 :  Normalized  power  spectral  density  of  signal  plus  colored  noise  (white  Gaussian 

noise  variance=10). 
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The  result  achieved  with  the  noise  matrix  eigenvector  transformation  is  shown  in 
Figure  4.8: 


Figure  4.8:  Normalized  power  spectral  density  of  whitened  signal  using 

eigendecomposition. 


Figure  4.8  shows  that  the  transformation  has  whitened  both  the  signal  and  the  noise 
components.  None  of  the  signal  remains  despite  low  noise  power,  and  detection  would  not 
be  possible.  The  reason  for  this  is  due  to  the  singular  condition  of  the  noise  autocorrelation 
matrix.  Although  the  rank  of  this  matrix  equals  256,  its  determinant  equals  zero.  This  result 
is  unexpected  and  may  be  a  function  of  the  magmtudes  of  the  values  in  the  correlation 
matrix.  The  decomposition  of  the  matrix  results  in  eigenvalues  which  are  less  than  one  and 

where  many  are  on  the  order  of  lO"®.  The  eigenvector  decomposition  method  is  particularly 
susceptible  to  this  problem.  The  singular  value  decomposition  technique  discussed  next 
provides  a  way  to  accurately  compute  the  eigenvectors  of  a  singular  matrix. 

4.  Singular  Value  Decomposition  (SVD) 

As  stated  earlier,  the  SVD  techmque  provides  an  alternative  way  to  compute  the 
matrix  of  eigenvectors  for  a  given  correlation  matrix  without  having  to  compute  it  first.  It 
does  so  by  using  a  matrix  of  data.  In  this  expaiment  the  data  used  is  that  of  the  colored 
noise.  It  is  assumed  that  sequences  of  recorded  noise  have  been  obtained  for  the  area  in 
question  and  that  this  noise  is  representative  of  the  actual  conditions.  The  noise  data  matrix 
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used  is  of  dimension  48  X  1024,  and  is  created  by  using  noise  data  of  different  realizations 
than  that  used  in  the  input  sequence  x.  From  this  noise  data  matrix,  the  SVD  algorithm 
computes  a  matrix  of  eigenvectors  and  a  matrix  of  singular  values.  The  matrix  of 
eigenvectors  is  then  used  in  the  same  way  as  for  the  previous  experiment  to  decorrelate  the 
noise  component  of  the  received  sequence.  The  figures  shown  were  computed  using  the 
Matlab  program  given  in  Appendix  K. 

Figure  4.9  shows  the  normalized  averaged  power  spectral  density  for  the  input 
sequence  whitened  via  the  SVD  decomposition.  In  this  case  the  noise  power  is  set  by  using 
a  white  Gaussian  noise  variance  of  50. 


As  can  be  observed,  the  noise  power  is  such  that  the  sinusoids  are  even  difficult  to 
distinguish  with  the  averaged  periodogram.  Even  the  use  of  a  bandpass  filters  around  the 
frequencies  of  interest  or  a  low  pass  filter  removing  the  stronger  noise  at  lower  frequencies 
would  make  detection  very  difficult. 

Figure  4.10  shows  the  averaged  power  spectral  density  of  the  transformed  sequence. 
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It  can  be  observed  that  the  results  are  far  better  with  the  SVD  method  than  with  the 
eigenvector  factorization  method  despite  the  fact  that  both  methods,  in  theory,  should 
arrive  at  the  matrix  of  eigenvector  of  the  same  correlation  matrix.  Sinusoid  #1  and  Sinusoid 
#2  are  readily  detectable.  Sinusoid  #3  however  only  appears  as  a  small  protrusion  and 
would  not  be  detectable.  This  result  demonstrates  the  superiority  of  the  SVD  algorithm 
when  the  correlation  matrix  of  the  noise  is  poorly  conditioned  or  singular.  As  stated  in 
Chapter  m,  the  SVD  method  is  far  more  robust  than  the  normal  eigenvector/eigenvalue 
decomposition  method  and  does  provide  the  means  to  achieve  optimal  detection  as  a 
prewhitener. 

5.  Mahalanobis  Transformation 

The  Mahalanobis  transformation  is  another  transformation  method  based  on  the 
eigendecomposition  of  the  theoretical  noise  correlation  matrix.  The  results  achieved  for 
this  transformation  are  shown  here.  These  results  were  computed  using  the  Matlab  program 
shown  in  Appendix  L.  The  results  where  achieved  by  separating  the  test  sequence  x  in  four 
blocks  and  transforming  each  of  these  according  to  Equation  3.96. 

Figure  4. 1 1  shows  the  normalized  power  spectral  density  for  the  input  sequence.  The 
noise  power  is  reduced  by  setting  the  white  Gaussian  noise  variance  to  25. 
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Figure  4. 1 1:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=25). 


Figure  4.12  shows  the  power  spectral  density  of  the  transformed  sequence. 
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Figure  4.12:  Normalized  power  spectral  density  of  whitened  signal  using  Mahalanobis 

transformation. 

It  can  be  observed  that  although  the  noise  has  maintained  some  correlation,  it  has  been 
considerably  whitened  when  compared  to  Figure  4.11.  At  the  higher  frequencies  where 
Sinusoid  #3  is  located,  the  transformation  has  maintamed  its  local  signal  to  noise  ratio  in  a 
way  in  which  it  would  be  easily  detected.  At  lower  frequencies  where  Sinusoid  #1  and  #2 
are,  the  signal  to  noise  ratio  at  that  location  along  the  frequency  axis  appears  to  have  been 
maintained  but  since  the  noise  floor  continues  rising  until  about  0.32 /s,  and  is  more  than 
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twice  the  height  at  that  frequency  than  it  is  at  0.15  fs,  detection  would  require  additinnal 
processing. 

6.  Lower-Upper  Factorization  (LDU) 

The  LDU  factorization  is  the  first  of  the  whitening  transformations  to  use  a  triangular 
matrix,  that  is  shown  here.  Figure  4.13  shows  the  normalized  power  spectral  density  for  the 
input  sequence.  In  this  case  the  noise  power  is  achieved  by  using  a  white  Gaussian  noise 
variance  of  25.  The  figures  shown  were  computed  using  the  Matlab  program  given  in 
Appendix  M. 


Figure  4.13:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=25). 

This  sequence  is  then  whitened  with  the  L  triangular  matrix  of  the  LDU  decomposition 
of  the  theoretical  correlation  matrix,  in  accordance  with  Equation  3.102. 

Figure  4.14  presents  the  averaged  power  spectral  density  of  the  whitened  sequence. 
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Figure  4.14:  Normalized  power  spectral  density  of  whitened  sequence  using  LDU. 


It  can  be  observed  that  the  transformation  has  not  succeeded  in  completely  whitening 
the  input  sequence.  From  0  to  0.02  there  is  an  upwards  curve  followed  by  a  steady 
descending  slope  from  0.025  to  0.5  U  However,  the  three  sinusoids  are  clearly 
distinguishable  with  few  occasional  peaks  which  could  lead  to  false  alarms.  Overall  the 
results  show  that  the  local  signal  to  noise  ratio  has  been  maintained  or  somewhat  improved 
but  that  the  noise  is  not  fully  whitened  especially  at  very  low  frequencies  where  the  slope 
of  the  noise  is  greatest. 

7.  Upper-Lower  Factorization  (UDL) 

This  transformation  is  essentially  the  same  as  above  but  it  uses  a  upper  triangular 
matrix.  This  leads  to  the  anti-causal  transformation  of  the  theoretical  noise  correlation 
matrix  shown  at  Equation  3.111.  The  results  of  this  whitening  on  the  same  input  sequence 
used  for  the  UDU  factorization,  shown  at  Figure  4.13,  is  shown  at  Figure  4.15. 


78 


1 


0  0.05  0.1  0.15  0.2  0.25  0.3  0.35  0.4  0.45  0.5 

Fraction  of  sampling  frequency 

Figure  4. 15:  Normalized  power  spectral  density  of  whitened  sequence  using  UDL. 

It  can  be  observed  that  the  results  achieved  with  the  UDL  transformation  is  very 
similar  to  that  achieved  with  the  LDU  transformation.  The  anti-causality  of  the 
transformation  seems  to  have  no  effect  on  the  performance  of  the  whitener. 

8.  Cholesky  Factorization 

As  stated  in  Chapter  DI,  the  Cholesky  decomposition  is  based  on  the  LDU  and  UDL 
decompositions.  As  such,  results  achieved  with  it  are  expected  to  be  similar  to  those 
resulting  from  those  methods.  As  is  done  for  the  LDU  and  UDL  decompositions,  a  white 
Gaussian  noise  variance  of  25  is  used  to  create  the  colored  noise  and  results  in  the  power 
spectral  density  for  the  input  sequence  shown  at  Figure  4. 16. 
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Figure  4.16:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=25). 

The  next  two  figures  present  m  turn  the  causal  and  anti-causal  Cholesky 
transformation  results,  related  respectively  to  the  LDU  and  UDL  factorizations.  Figure 
4.17,  shows  the  power  spectral  density  of  the  transformed  sequence  using  the  Cholesky 
factor  consisting  of  a  lower  triangular  matrix,  therefore  resulting  in  a  causal  transformation. 
This  transformation  is  performed  using  Equation  3.115.  The  Matlab  program  used  to 
perform  the  transformations  and  generate  the  figures  is  included  at  Appendix  N. 
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Figure  4.17:  Normalized  power  spectral  density  of  whitened  sequence  using  lower 

triangular  Cholesky  factor. 


Figure  4.18  presents  the  averaged  power  spectral  density  of  the  transformed  sequence 
using  the  upper  triangular  Cholesky  matrix.  This  transformation  is  anti-causal. 
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We  can  see  that  the  results  achieved  are  quite  good  in  that  the  noise  is  whitened  and 
the  ratio  of  the  spectral  heights  of  each  of  the  sinusoids  and  the  noise  at  the  respective 
spectral  locations  are  also  maintained.  Some  spurious  peaks  may  cause  false  alarms  for 
Sinusoid  #2.  The  causaMty  or  anti-causality  of  the  transformation  seems  to  have  no  effect 
on  the  performance  achieved.  The  three  sinusoids  are  reliably  detected  and  have  maintained 
or  improved  their  local  SNR  relative  to  the  input  sequence.  Since  the  Cholesky 
decomposition  differs  from  the  LDU  decomposition  in  that  the  square  root  of  the  ftiagnnfti 
matrix  Dj^  scales  the  matrix  L ,  this  last  factor  appears  to  allow  the  method  to  correctly 
whiten  the  noise. 

9.  QR  Factorization 

QR  factorization,  like  SVD  makes  use  of  a  matrix  of  data  to  arrive  at  the  lower 
triangular  matrix  L  of  the  LDU  factorization.  In  this  case,  a  square  noise  only  data  matrix 
of  dimensions  256  X  256  is  used.  In  the  same  way  as  for  the  SVD  method,  the  noise  used 
in  the  data  matrix  is  independent  of  that  used  in  the  sequence  x.  Four  block  transforms  of 
length  256  are  joined  to  form  the  whitened  sequence  of  1024  points.  The  Matlab  program 
used  to  generate  the  results  is  shown  at  Appendix  O. 
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Figure  4.19  shows  the  power  spectral  density  for  the  input  sequence.  The  white 
Gaussian  noise  variance  is  set  to  50. 
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Figure  4.19:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=50. 


One  point  to  note  concerning  Figure  4.19  is  that  Sinusoid  #2  is  not  well  resolved  in 
this  realization  of  the  process.  This  result  is  also  observed  in  the  whitened  output  shown  in 
Figure  4.20. 
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Figure  4.20:  Normalized  power  spectral  density  of  whitened  signal  using  the  QR 
factorization  to  achieve  an  estimate  of  the  L  matrice.  (256X256) 

Figure  4.20  shows  that  the  colored  noise  is  successfully  whitened.  AH  three  sinusoids 
are  identifiable  although  the  variance  of  the  estimate  may  still  lead  to  some  false  alarms 
depending  on  the  height  of  the  detection  threshold.  As  stated  above,  the  Sinusoid  #2  peak 
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is  not  as  detectable  as  Sinusoid  #1  and  Sinusoid  #3  but  this  is  also  the  case  with  the  input 
sequence.  The  local  signal  to  noise  ratio  of  all  three  sinusoids  has  also  been  improved 
compared  to  the  averaged  periodogram  shown  in  Figure  4. 19.  This  result  seems  to  indicate 
that  the  numerical  accuracy  required  to  perform  successful  whitening  with  the 
transformations  discussed  herein  is  important  and  is  not  achievable  from  the  use  of  the 
theoretical  noise  correlation  matrix.  Overall  this  is  the  most  successful  whitener  among 
those  using  matrix  decomposition  methods  in  terms  of  truly  whitening  the  colored  noise. 

C.  RESULTS  OF  OTHER  NON-OPTIMUM  FILTERS 

1.  Differential  Filter 

This  section  shows  the  results  of  using  the  simple  differentiation  function  to  attempt 
whitening  as  discussed  in  Chapter  HI.  The  advantage  in  using  this  method,  should  it  be 
successful,  lies  in  its  simplicity  and  efficiency.  It  does  not  require  any  matrix  factorization 
process  nor  does  it  demand  numerous  multiplications,  but  it  can  only  be  used  in  certain 

2 

colored  noise  situations  (i.e.,  1// ,  1//  noise,  etc.).  The  figures  shown  were  computed 
using  the  Matlab  program  given  in  Appendix  H. 

The  averaged  periodogram  for  the  input  sequence  of  signal  plus  noise  is  shown  in 
Figure  4.21.  For  this  experiment  the  white  Gaussian  noise  variance  is  set  to  25. 


Fraction  of  sampling  frequency 


Figure  4.21:  Normalized  power  spectral  density  of  test  sequence  (white  Gaussian  noise 

variance=25). 
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Three  figures  are  shown  next.  They,  in  order,  show  the  first,  second  and  third  order 
differentiation  in  time  of  the  input  sequence.  As  discussed  in  Chapter  HI,  it  is  expected  that 
the  first  order  differentiation  of  the  sequence  should  provide  the  best  prewhitening.  Figure 
4.22  shows  the  averaged  periodogram  of  the  first  order  differentiation  of  the  input 
sequence. 
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Figure  4.22:  Normalized  power  spectral  density  of  whitened  sequence  using  first  order 

differentiation. 

As  with  some  of  the  other  whitening  techniques,  it  can  be  readily  observed  that 
although  the  noise  has  been  equalized  somewhat,  it  still  remains  colored.  Nevertheless,  we 
can  see  that  the  three  signal  peaks  are  distinct  and  detectable.  However,  the  form  of  the 
noise  power  spectral  density  would  lead  to  better  detection  results  in  using  a  white  noise 
approximation  between  the  frequencies  of  0.075 fs  -  0.2  fs.  In  that  area,  the  noise  spectrum 
is  relatively  flat  and  signals  such  as  Sinusoid  #2  produce  sharp  easily  detected  peaks. 
Although  Sinusoid  #1  and  Sinusoid  #3  are  not  in  areas  where  the  noise  power  has  leveled 
off,  the  local  signal  to  noise  ratio  of  all  three  sinusoids  appears  to  be  at  least  as  high  as  the 
local  signal  to  noise  ratio  shown  on  Figure  4.21.  All  three  sinusoids  of  Figure  4.22  account 
for  the  global  peaks  whereas  on  Figure  4.21,  Sinusoid  #2  and  Sinusoid  #3  have  less  powca- 
than  die  noise  in  the  low-pass  region. 

Because  the  noise  is  not  completely  equalized  by  using  the  first  order  differentiation, 
a  second  and  third  experiment  are  conducted  to  see  what  the  effect  of  second  order  and  third 
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order  differentiation  would  have  on  the  input  sequence.  Figure  4.23  shows  the  resulting 
averaged  periodogram  for  the  second  order  differentiation  and  Figure  4.24  contains  the 
results  of  the  third  order  differentiation. 


Figure  4.23:  Normalized  power  spectral  density  of  whitened  sequence  using  second  order 

differentiation. 


Figure  4.24:  Normalized  power  spectral  density  of  whitened  signal  using  third  order 

differentiation. 

The  effect  of  the  additional  orders  of  differentiation  can  be  readily  seen  on  the  two 
figures.  In  all  cases  the  noise  power  forms  a  characteristic  hump  which  begins  and  ends 
with  magnitudes  close  to  or  equal  to  zero.  For  the  first  order  differentiation,  the  top  of  the 
hump  is  centered  at  a  frequency  of  about  0.12/5.  For  the  second  order  differentiation,  the 
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hump  moved  to  higher  frequencies  and  its  summit  is  around  0.26  fs.  In  the  case  shown  in 
Figure  4.24,  the  somewhat  flat  noise  power  area  is  in  the  frequency  range  of  0.2/j-  0,32/s 
and  Sinusoid  #3  at  0.3  fs  is  clearly  detectable.  Sinusoid  #1  at  0.05  fs  is  far  lower  than  the 
noise  present  at  most  frequencies  and  is  not  detectable.  Sinusoid  #2  is  in  the  middle  of  the 
slope  leading  to  the  maximum  noise  power  and  has  maintained  its  local  signal  to  noise  ratio 
but  the  noise  at  higher  frequencies  has  higher  power  and  would  need  to  be  filtered  to  allow 
for  correct  detection.  For  the  third  order  differentiation,  the  hump  moves  to  higher 
frequency  with  the  top  of  the  noise  hump  situated  at  a  frequency  around  0.32  fs.  Sinusoid 
#2  and  Sinusoid  #3  remain  identifiable,  but  Sinusoid  #1  has  disappeared  along  with  the 
noise  in  the  frequencies  around  0.05  fs  and  lower.  Because  of  its  location  in  the  flat  portion 
of  the  noise  spectrum  Sinusoid  #3  again  is  sharply  defined  and  his  easily  detectable  witfi  a 
white  noise  assumption. 

This  method  is  a  very  easy  and  inexpensive  way  to  minimize  the  effects  of  low 
frequency  noise  that  has  spectral  roll-off  as  is  found  in  the  ocean  environment. 

2.  Prediction-Error  Filter  (PEF) 

This  section  demonstrates  the  whitening  property  of  the  prediction-error  filter.  As 
discussed  in  Chapter  HI,  the  success  of  the  prediction  error  filter  is  predicated  upon  the 
correlation  feature  of  the  signal  and  of  the  noise.  The  signal  used  for  this  experiment  is 
different  than  that  used  for  the  other  experiments  and  is: 

I] 

The  signal  is  changed  to  better  illustrate  the  performance  of  the  filter  and  its 
dependence  upon  the  correlation  properties  of  the  signal  and  colored  noise.  Figures  4.25 
and  4.26  show  the  theoretical  correlation  functions  of  the  colored  noise  and  of  the  signal 
for  the  first  21  lags. 


5[/i]  =  0.475cos|^27in|  +  0.15cos^2nn|j  +  0.1cos^27i/i 
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Figure  4.25:  Autocorrelation  function  of  the  Butterworth  low-pass  noise:  First  21  points. 

Figure  4.25  shows  that  the  magnitude  of  the  autocorrelation  function  of  the  noise 
rapidly  decreases  with  the  number  of  delays.  The  autocorrelation  function  of  underwater 
noise  is  very  similar  because  it  also  has  high  noise  power  for  low  frequencies.  This  means 
that  if  a  delay  must  be  chosen  to  maximize  the  correlation  of  the  noise,  then  this  delay 
should  be  small.  If  the  noise  is  white,  the  magnitude  would  be  non-zero  only  at  a  delay  of 
zero. 

Figure  4.26  shows  the  autocorrelation  function  of  the  signal. 


Figure  4.26:  Autocorrelation  function  of  the  signal:  s=.475*cos(2*pi*t/8) 
-K0.15*cos(2*pi*t*2/8)  +  0.1*cos(2*pi*t*3/8) 


As  expected,  the  autocorrelation  function  of  the  sinusoidal  signal  is  also  sinusoidal. 
We  note  that  the  signal  is  formed  from  three  sinusoids  each  of  which  has  its  own  correlation 
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function.  The  net  results  shown  here  presents  the  sum  of  all  three  autocorrelation  functions. 
Therefore,  if  a  delay  is  chosen  where  the  autocorrelation  function  of  the  signal  does  not 
have  a  magmtude  of  zero  as  shown  on  Figure  4.26,  each  of  the  three  sinusoids  is  affected 
differently  based  on  the  magnitude  of  their  respective  autocorrelation  function  at  that  delay. 

In  order  to  demonstrate  the  whitening  properties  of  the  prediction-error  filter,  an 
adaptive  order  32  LMS  filter  was  constructed.  Data  is  input  to  the  filter  until  the  weights 
have  converged  and  all  results  are  computed  only  with  data  obtained  after  convergence. 
The  whitened  sequence  is  the  error  between  the  predicted  value  out  of  the  predictor,  and 
the  actual  value  of  the  sequence.  Therefore,  if  the  signals  are  to  be  detected  in  the  error,  the 
delay  d  (see  Figure  3.9)  must  be  chosen  such  that  the  correlation  of  the  noise  is  maximized 
while  the  correlation  of  the  signal  is  minimized. 

The  next  three  figures  present  the  results  of  the  experiment.  The  Matlab  code  of 
Appendix  Q  was  used  to  compute  these  results. 

Figure  4.27  shows  an  averaged  periodogram  of  the  input  sequence.  The  white 
Gaussian  noise  variance  used  to  create  the  colored  noise  is  100.  Experimentation  shows 
that  the  correlation  of  the  noise  samples  is  most  sensitive  to  delays.  A  delay  of  one  achieved 
the  best  whitening.  We  expect  that  this  is  due  to  the  fact  that  the  colored  noise  has  a 
correlation  function  which  asymptotically  decreases  with  delay.  Therefore,  as  the  delay  is 
increased,  the  predictor  quickly  loses  its  ability  to  estimate  the  noise  and  therefore  remove 
it  from  the  error  sequence.  The  signals  can  be  seen  0.125/*,  0.25/*,  and,  0.375/*. 
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Fraction  of  sampling  frequency 


Figure  4.27:  Normalized  power  spectral  density  of  signal  plus  noise  (white  Gaussian  noise 

variance=l(X)). 


Figure  4.28  shows  the  power  spectral  density  of  the  output  of  the  predictor. 


Fraction  of  sampling  frequency 

Figure  4.28:  Normalized  power  spectral  density  of  estimated  signal 

For  a  delay  of  one,  the  autocorrelation  of  none  of  the  three  sinusoids  is  zero.  Because 
of  this,  the  predictor  is  also  able  to  estimate  their  presence  as  well  as  that  of  the  colored 
noise.  Since  the  delay  chosen  is  short,  it  can  be  expected  that  low  frequency  noise 
components  and  signals  at  low  frequencies  would  be  better  predicted  and  tho'efore  more 
easily  removed  from  the  whitened  sequence  or  the  error  sequence.  In  the  case  of  signals  in 
low-pass  noise,  this  is  a  good  feature  since  the  noise  at  low  frequencies  has  far  greater 
power  than  the  signals  which  are  at  higher  spectral  locations.  We  note  that  the  averaged 
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petiodogram  of  the  predicted  signal  looks  very  similar  to  the  averaged  periodogram  of  the 
input  sequence. 

Figure  4.29  shows  the  averaged  periodogram  of  the  error  containing  the  whitened 
sequence. 


Figure  4.29:  Normalized  power  spectral  density  of  predicted  error:  whitened  signal. 


It  can  be  seen  that  the  noise  has  been  considerably  whitened  by  the  PEF.  At 
frequencies  below  0.27/„  the  noise  is  effectively  whitened.  At  higher  frequencies,  the  small 
delay  appears  to  have  partly  decorrelated  some  of  the  noise  power,  and  some  coloration 
remains.  The  three  signals  also  appear  clearly  despite  a  delay  that  did  not  completely 
decorrelate  them  as  indicated  by  their  presence  in  the  predicted  sequence.  Comparing 
Figures  4.27  and  4.29  shows  that  the  ratios  between  the  spectral  heights  of  each  sinusoid 
and  the  noise  at  their  respective  spectral  location  have  been  preserved.  Different  delays 
however,  have  a  serious  effect  of  the  performance  of  this  whitener.  As  the  selected  delay 
increases,  the  noise  becomes  more  decorrelated,  and  the  predictor  is  less  able  to  produce 
correct  estimates.  This  then  results  in  more  of  the  colored  noise  to  be  part  of  the  error. 

When  this  technique  is  used  for  detection  in  colored  noise,  its  performance  is 
dependent  upon  the  correlation  properties  of  the  signal  and  noise.  In  the  case  of  the  noise, 
imderwater  ambient  noise  is  always  best  predicted  by  choosing  small  delays.  For  the  signal, 
which  must  be  decorrelated  as  much  as  possible,  a  study  of  its  correlation  properties  must 
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be  undertaken  to  ensure  that  the  delay  chosen  does  not  occur  at  a  point  where  the  signal 
autocorrelation  function  is  at  its  maximum.  Should  this  occur,  the  predictor  is  then 
completely  able  to  predict  the  signal  in  its  estimate  and  therefore  none  of  the  signal  appears 
in  the  error  sequence. 

An  advantage  of  the  method  presented  herein,  is  that  it  is  adaptive.  As  sea  going 
vessels  are  underway,  the  noise  characteristic  of  the  imderwater  environment  changes.  It  is 
reasonable  to  assume  that  this  could  affect  our  choice  of  delay.  However,  the  variations  in 
the  ambient  noise  spectra  usually  only  minimally  affect  this  consideration,  because  the 
primary  effect  resulting  in  much  higher  noise  levels  at  low  frequency  is  the  sound 
attenuation  in  sea  water.  Therefore  using  an  adaptive  filter  in  this  manner  allows  these 
changes  to  be  considered  in  the  detection  process  without  modifying  the  experimental 
parameters. 
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V.  CONCLUSIONS 


In  this  thesis,  we  have  presented  a  number  of  techniques  with  which  to  perform 
discrete  time  detection  of  signals  in  colored  noise.  Ute  matched  filter  for  colored  noise  is 
derived  and  is  shown  to  produce  excellent  results.  Whrai  the  correlation  matrix  of  the  noise 
can  be  estimated,  the  matched  filter  for  colored  noise  provides  the  optimal  detection 
solution.  The  derivation  for  the  Maximum  A  Postaiori  oiteiion  for  colored  noise  is  also 
shown.  The  inverse  filter  is  demonstrated  to  be  a  very  efficient  way  to  whiten  the  noise  by 
effectively  equalizing  the  colored  noise  by  the  use  of  its  inverse  power  spectral  density. 
This  method  is  simple  and  effective  with  the  only  requiimnent  being  that  the  power  spectral 
density  of  the  noise  be  known.  Several  matrix  decmnposition  techniques  are  addressed. 
These  techniques  are  based  on  the  use  of  a  block  transformation  ^proach  to  whitening  a 
sequence.  They  generally  use  a  matrix  of  orthonmmal  vectors  derived  from  the  noise,  to 
transform  a  data  sequence  in  order  to  decorrelate  die  noise  component  of  that  sequence.  Tf 
the  colored  noise  can  be  decorrelated,  then  detection  in  white  noise  can  be  performed.  The 
procedures  to  obtain  a  matrix  of  orthonormal  vecUws  rely  on  a  factcnization  of  either  the 
correlation,  or  covariance  matrix  of  the  colored  noise  process,  ot  of  a  noise  data  matrix.  The 
methods  using  the  data  matrix  (i.e.,  signal  free  noise  recording)  are  less  susceptible  to  a 
poorly  conditioned  or  singular  correlation  matrix  whidi  is  necessarily  estimated.  They  also 
result  in  better  numerical  accuracy  which  was  found  to  be  important  for  the  whitening 
process.  Some  of  the  methods  examined,  for  example,  flic  LDU,  DDL,  Cholesky,  and,  QR 
factorizations,  form  a  matrix  which  is  triangular.  Alfliough  some  matrix  decomposition 
methods  require  less  processing  resources  than  others,  all  of  them  remain  vary  inefficient 
compared  to  other  available  techniques  such  as  the  matched  filter  and  the  inverse  filter. 
Some  of  these  techniques,  although  promising  in  theory,  are  sometimes  unable  to 
completely  whiten  the  colored  noise  as  predicted.  For  flie  low-pass  Butterworth  colored 
noise  used,  the  eigenvector  factorization  method  was  unsuccessful  in  allowing  for  the 
detection  of  any  of  the  three  sinusoids.  As  a  resulte  of  flie  poor  conditioning  of  the  noise 
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autocorrelation  matrix,  the  transformation  using  the  matrix  of  eigenvectors  also  whitens  the 
signal  along  with  the  noise.  Three  methods  work  well  in  preserving  the  signal  to  noise  ratio 
of  the  sinusoids  while  whitening  the  noise.  The  singular  value  decomposition  computes  a 
more  accurate  matrix  of  eigenvectors  for  the  noise  autocorrelation  matrix  and  functions 
well  in  spite  of  poor  conditioning  of  the  autocorrelation  matrix.  It  whitens  the  noise  and 
generally  maintains  signal  power.  The  QR  technique  also  uses  a  noise  data  matrix  to 
achieve  a  more  precise  estimate  of  the  L  matrix  of  the  LDU  factorization  of  the  correlation 
matrix.  \s  shown  in  Chapter  4,  the  QR  method  effectively  decorrelates  the  coined  noise 
and  allows  the  use  of  an  optimal  white  noise  detector,  to  achieve  optimal  detection.  The 
Cholesky  factorization  is  able  to  correctly  whiten  the  noise  and  maintflin  the  appropriate 
spectral  heights  ratio  but  is  not  as  capable  as  the  QR  factorization  when  the  input  noise 
power  is  greater.  The  differential  operator  can  be  a  very  efficient  whitener  when  the  colored 
noise  is  inversely  proportional  to  a  power  of  the  frequency.  The  whitening  propaty  of  the 
prediction-error  filter  is  demonstrated  and  is  shown  to  be  effective.  However,  it  is  strongly 
dependent  on  the  correlation  properties  of  the  signals  to  be  detected.  Some  periodic  signals 
may  not  be  detected  if  the  chosen  delay  occurs  at  a  maximum  of  the  signal  correlation 
function.  When  the  delay  is  at  a  point  of  high  correlation  for  the  noise,  and  at  a  point  of  low 
correlation  for  the  signal,  the  prediction  error  filter  can  effectively  whiten  the  noise  and 
provide  good  ouq>ut  signal  to  noise  ratios. 
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APPENDIX  A  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  22,  “THERMAL  AGITATION  SOUND  PRESSURE  LEVEL 

(SPL)”. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Plots  thermal  agitation  noise  levels 
% 

%  Author:  Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  tiiermalnoisem 


f=logspace(0,7,500);%frequency  range 

SPL=-101+20*logl0(f); 

semilogx(f,SPL) 
axis([l  10^7-100  40]) 
grid 

xlabel( ‘Frequency  (Hz)’) 

ylabel(‘SPL  (dB)  re  2X10^3  uPa  for  a  1  Hz  Bandwidth’) 
print  thermalBg 

!ps3epsi  <theimalBg.ps  >thermalfig.epsi 
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APPENDIX  B  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  2.4,  “COMPARISON  OF  ABSORPTION  COEFFICIENTS  IN 
SEA  WATER,  DUE  TO  SHEAR  AND  VOLUME  VISCOSITY,  AND 

SHEAR  VISCOSITY  ALONE”. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Thesis:  Plot  of  viscosity  absorption  coefficients  iaw  Fisher  and  Simmons 
%salinity:  35  ppt 
%PH=8.0 

%ref  page  104  Urick 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  absorptionf21.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clf 

Mogspace(0,7,500);%frequency  range 
T=5;  %temperature  in  Celcius 

Po=l;  %pressure  in  atm 

TK=T+273;%temp  in  Kelvin 

fl=1.32*10^3*TK*exp(-1700/TK);  %relaxation  freq  of  boric  add 
f2=1.55*10^7*TK*exp(-3052/TK);%relaxation  freq  of  MgS04 

%  for  temps  and  pressures  in  ranges  0-30iC  and  1-400  atm 
A=8.95*10^(-8)*(1+2.3*10^(-2)*T-5.1*10^(-4)*T'^2); 

B=4.88*  l(y'(-7)*(  1+1 .3*  10^(-2)*T)*(  1-0.9*  10^-3)*Po); 
C=4.76*10^(-13)*(l-4*10^(-2)*T+5.9*10A(-4)*T'^2)*(l-3.8*10^(-4)*Po); 

%attenuation  coeff  in  dB/m 
attboric=A*f  1  *f.''2./(f  l''2+f .'^2); 
attMgS04=B*f2*f.A2./(f2A2+f.A2); 
attwater=C*f.'^2; 

atb=attboric+attMgS04+attwater; 

mus=.01;%shear  viscosity 
muv=2.81*mus;%volume  viscosity 
density=  1  ;%gm/cm3 
c=1.5*10^5;%cm/s 
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%pure  water 

attshearvisc=16*pi''2*mus*f.''2./(3*density*c''3); 

attvisc=16*pi^2*f.'^2*(mus+0.75*muv)/(3*density*c^3); 

loglog(f,att,’-’,f,1000*attshearvisc,’-’4,attwater,’-.’) 
axis([l  10^7  lO^(-lO)  10^2]) 
ylabel(‘Absorption  coefficient  (dB/m)’) 
xlabel(‘Frequency  (Hz)’) 

h=legend(‘  sea  water’,’  shear  and  volume’,’  shear  (theoretical)’) 
axes(h) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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APPENDIX  C  ■  MATLAB  PROGRAM  USED  TO  COMPUTE  FIGURE 
2.5,  “SOUND  ABSORPTION  AT  5°C,  1  ATM,  35  PPT  SALINITY,  AND 
PH=8.0  (ACCORDING  TO  FISHER  AND  SIMMONS)”. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Thesis:  Plot  of  absorption  coefficients  according  to  Fisher  and  Simmons 
%salinity:  35  ppt 
%PH=8.0 

%ref  page  158  KFCS 
%  For  second  figure  Chap  2 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  absorptionf22.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clf 

clear 

f=logspace(0,7,500);%ffequency  range 
T=5;%temperature  in  Celcius 
Po=l;%pressure  in  atm 

TK=T+273;%temp  in  Kelvin 

fl=1.32*10''3*TK*exp(-17(X)/rK);  %relaxation  freq  of  boric  acid 
f2=1.55*10^7*TK*exp(-3052/rK);%relaxation  ffeq  of  MgS04 

%  for  temps  and  pressures  in  ranges  0-30jC  and  1-400  atm 

A=8.95*10^(-8)*(1+2.3*10^(-2)*T-5.1*10^(-4)*T''2); 

B=4.88*10'^(-7)*(l+1.3*10A(-2)*T)*(l-0.9*l(yK-3)*Po); 

C=4.76*10^(-13)*(l-4*10^(-2)*T+5.9*10^(-4)*T'^2)*(l-3.8*10^(-4)*Po); 

%attenuation  coeff  in  dB/m 

attboric=A*f  1  ♦f A2./(f  1'^2+f  ^2); 

attMgS04=B*f2*fA2./(f2A2+fA2); 

attwater=C*fA2; 

att=attboric+attM^S04+attwater; 
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loglog(f,att,’-’,f,attwater,’--’Xattboric,’:’,f,attMgS04,’-.’) 
axis([l  10^7  Ky^C-lO)  10^2]) 
ylabel(‘Absorption  coefficient  (dB/m)’) 
xlabel(‘Frequency  (Hz)’) 

h=legend(‘  seawater’,’  fresh  water’,’  boric  acid’,’  magnesium  sulfate’) 
axes(h) 


APPENDIX  D  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  2.7,  “NORMALIZED  AVERAGE  SOUND  PRESSURE 
LEVEL  OF  AMBIENT  NOISE  AT  POINT  SUR  SOSUS  ARRAY 
(SAMPLING  FREQUENCY  8  KHZ)” 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%PSD  of  ocean  noise  from  sumoise  file  taken  at  8  kHz  on  a  log  scale 
%  From  sumoise.au 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  noisesurPSDlogpad.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clf 

orient  tall 

acnoise=auread(‘sumoise.au’); 

%nacnoise=acnoise-mean(acnoise);  %remove  DC 
index=l; 

for  1=1:2048:40960 

PSD(:,index)=20*logl0(abs(fft(acnoise(i:i+2047),8192))); 

index=dndex+l; 

end 

PSDave=sum(PSD’)./(index- 1); 

%plot  in  semilog  space 
clf 

subplot(211) 

point=logspace(0,logl0(4096),500); 
semilogx(point*8000./8192,PSDave(point)-max(PSDave(point))) 
xlabel( ‘Frequency  (Hz)’) 
ylabel(‘Magnitude  (dB)’) 

axis([10'^l  10Mmin(PSDave(point)-max(PSDave(point)))  0]) 
print  fsumoisePSD 

IpsSepsi  <fsumoisePSD.ps  >fsumoise.epsi 
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APPENDIX  E  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  3.4,  “DETECTION  REGIONS:  GAUSSIAN  PROBABILITY 

DENSITY  FUNCTIONS”. 


%  Program  that  plots  detection  regions 
% 


%  Author:Capt  M.A.  Cloutier 
%Date:  26  Oct  95 
%  Name:  gauss.m 


orient  taU 


i=-15:.l:25; 

m=10; 

var=4; 

pl=l/(sqrt(2*pi*var'^2))*exp(-(i-m).''2./(2*var'^2)); 

pO=l/(sqrt(2*pi*var^2))*exp(-(i)A2./(2*var'^2)); 


T=linspace(0,0.75/(sqrt(2*pi*var'^2)),20); 

j=5*ones(si2e(T)); 


subplot(311) 
plot(i,pU,pOo,T) 
hold  on 


fill([i(202)  i(202:401)],[0  pl(202:401)].’m’); 
fill([i(202)  i(202:401)],[0  p0(202:401)]/g’); 

fm([i(l:200)  i(200)].[p0(l:200)  0].’c’); 
fill([i(l:200)  i(200)],[pl(l:200)  0].’y’); 

axisC'ofT) 
print  gaussfig 


!ps3epsi  <gaussfig.ps  >gaussfig.epsi 


APPENDIX  F  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  3.10,  “CORRELATION  FUNCTION  OF  FIRST  ORDER 
LOW-PASS  BUTTERWORTH  COLORED  NOISE”. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Program  to  compute  the  theoretical  autocorr  of  noise  from  the  PSD 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  2  Dec  95 
%  Name:  wienerk.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%7o% 

clf 

orient  tall 

1=1024; 

t=l:l; 

s=  l*sin(2*pi*t*l/4);%f  is  at  0.25  fs 

%varia=input(‘What  is  the  variance  of  the  white  noise?  ‘); 
varia=l; 

w=sqrt(varia)*randn(  1  ,length(s));%white  noise 

[B,A]=butter(l,0.1);%color  the  white  noise 
n=filter(B,A,w); 


x=s+n; 

PSD=abs(fft(n,2048)).'^2./2048;%PSD  of  the  noise 
[H,P]=freqz(B,A,2047);%H  contains  the  transfer  function 

Rnt=real(ifft(abs([H.’  fliplr(H.’)]).''2));%with  Wiener-Khinchine  theorem  find 
correlation  matrix 

Rn=xcorr(n);  %has  length  2*1-1 
Rnl=Rn.’;%form  row  vector  instead  of  column  vector 
%form  new  Rn  with  peak  at  zero  instead  of  centered 
Rn2=[Rnl(1024:2047)  Rnl(l:1023)]; 
impulse=dimpulse(B ,  A. 1024); 

Rnti=xco]T(impulse); 
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Rnmat=Rnti( 1024:2047); 
subplot(41 1) 

plot(0:1023,Rn2(l;1024)/1024,’-’,  0:1023, Rnt(l:1024),’.’) 
title(‘Plots  of  Rn  with  XCORR  vs  Rn  with  Wiener-Khinchine’) 

subplot(412) 

plot(0:19,Rn2(l:20)/1024, 0:19,  Rnt(l:20),0:19,Rnti(1024:1043)) 
title(‘Plots  of  first  20  points  of  correlations’) 

subplot(413) 

plot(abs(fft(Rn2/2047))) 

title(‘PSDwithRn2’) 

subplot(414) 

plot(abs(fft(Rnt))) 

title(‘PSD  with  Rnt’) 

print  figwk 

pause 

clf 

%for  thesis  chap4 
subplot(41 1) 

plot(0:25,Rn2(l:26)/1024,’-’.0:25,zeros(l,26),’.’) 

axis([025-.05.15]) 

xlabel(‘Delay’) 

ylabel(‘Amplitude’) 

print  figcorr 

!ps3epsi  <figcorr.ps  >figcorr.epsi 


APPENDIX  G  -  MATLAB  PROGRAM  USED  TO  GENERATE 
FIGURE  4.1,  “NORMALIZED  AVERAGE  POWER  SPECTRAL 
DENSITY  OF  SIGNAL  PLUS  NOISE”. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  This  will  create  the  signal  plus  noise  signal 
%  and  will  produce  the  diesis  figure  in  Chap  4 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  signalf.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

clear 

clf 

orient  landscape 
%  generate  data 

numb=100;%mmiba’  of  realizations  for  the  averaging 

l=numb*1024; 

t=l:l; 

s=  .475*cos(2*pi*t/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 

%varia=input(‘What  is  the  variance  of  the  white  noise?  ‘); 

varia=100; 

w=sqrt(varia)*randn(l,length(s));%white  noise 

[B,A]=butter(l,0.1);%color  the  white  noise 
colored=filter(B,A,w); 

x=s+colored; 

PSD=zeros(numb,2048); 

index=l; 

for  i=l:numb 

PSD(i,:)=abs(fft(x(index:(index+1023)),2048))A2./2048; 
index=index+ 1024; 
end 

PSDxave=sum(PSD)./numb; 
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PSDxave=PSDxave./max(PSDxave);%nonnalize 


%  Output  results 


plot(( 1 : 1024)/1024/2,PSDxave(l :  1024)); 
xlabel( ‘Fraction  of  sampling  frequency’) 
ylabel(‘Average(abs(fft(x,2048)).'^2)’) 

text(.08,.95*max(PSDxave),’Sin  #1’) 

axes(‘position’,[.4  .35  .15  .5]) 

plot((261:360)/1024/2,PSDxave(261:360)) 

title(‘Sin#2’) 

axis([261/2/1024  360/2/1024  0 .2]) 
grid 

axes(‘position’,[.7  .35  .15  .5]) 

plot((57 1 :670)/1024/2,PSDxave(57 1 :670)) 

title(‘Sin#3’) 

axis([57 1/2/1024  670/2/1024  0  .025]) 
grid 


print  sigffigl 

!ps3epsi  <sigffigl.ps  >sigffig.epsi 


APPENDIX  H  -  MATLAB  PROGRAM  IMPLEMENTING  THE 
DISCRETE  TIME  MATCHED  FILTER  FOR  COLORED  NOISE 


%  Discrete  Matched  Filter  for  colored  noise 


%  AuthonCapt  M.A.  Cloutier 
%  Date:  5  Dec  1995 
%  Name:  mffin.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

orient  tall 

%  generate  data 

numb=10;%number  of  realizations  for  the  averaging 

l=numb*1024; 

t=l:l; 

s=  .475*cos(2*pi*l/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 

%varia=dnput(‘What  is  the  variance  of  the  white  noise?  ‘)l 

varia=100; 

w=sqrt(varia)*randn(Uength(s));%white  noise 

[B  A]=butter(l,0.1);%color  the  white  noise 
colored=filter(B,A.w); 

x=s+colored; 

x=Teshape(x,  1024^umb)’; 
colored=reshape(colored, 1024,ninnb) 

%  Find  correlation  matrix  of  the  noise  from  PSD 

[H,P]=freqz(B,A2047);%H  contains  the  transfer  function 

Rn=real(ifft(abs([H.’  fliplr(H.’)])A2));%with  Wiener-Khinchine  theorem  find 
correlation  matrix 


Rntoep=toeplitz(Rn(  1 : 128)); 


%calculate  h 
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h=mv(Rntoep)*£liplr(s(  1 : 1 28)  ’)  ./sqrt(s(  1 : 128)*inv(Rntoep)*s(  1 : 1 28)  ’ ); 
y=zeros(numb, 1024); 

PSDy=zeros(numb, 1024) ; 

PSDx=zeros(numb,1024); 

fori=l;numb 

forj=l:(1024-127) 

y(ij)=h/*fliplr(x(i,j:a+127)))’; 

end 

PSDy(i.:)=(abs(fft(y(i.:).1024))A2); 

PSDx(i,:)=(abs(fft(x(i,:),1024))A2); 

end 

PSDyave=sum(PSDy)./numb; 

PSDyave=PSDyave/max(PSDyave); 

PSDxave=sum(PSDx)./nunib; 

PSDxave=PSDxave/niax(PSDxave); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%output  the  estimates 

clf 

subplot(311) 

plot((0:51 1)./1024,  PSDxave(l:512)) 
ylabel(‘Averaged(abs(fft(x))  A2)’) 
xlabel(‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot((130: 180)/1024,PSDxave(130: 180)./max(PSDxave(130: 180))) 
axis([130/1024  180/10240  1]) 

axes(‘position’,[.7  .78  .15  .12]) 

plot((285:335)/1024,PSDxave(285:335)/max(PSDxave(285:335))) 
axis([285/1024  335/10240  1]) 

print  mffrgx 

!ps3epsi  <mffigx.ps  >mffigx.epsi 
clf 

subplot(311) 
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plot((0:5 1 1)./1024,  PSDyave(l  :5 12)) 
ylabel(‘Averaged(abs(fft(y)).'^2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  mffigy 

!ps3epsi  <mffigy.ps  >mffigy.epsi 


APPENDIX  I  -  MATLAB  PROGRAM  IMPLEMENTING  THE 


INVERSE  FILTER 


%  Program  to  compute  coloring  and  whitening  of  a  sequence  with  white  noise 
%  using  an  inverse  filter 


% 


%  Author:Capt  M,A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  invfln.m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

orient  tall 


numb=10; 

P=zeros(numb,5 12); 

Pxa=zeros(numb4 12); 

for  ind=l:numb 

1=1024; 

t=l:l; 

s=.475*cos(2*pi*t/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10); 
varia=25; 

w=sqrt(varia)*randn(l,length(s));%white  noise 

[B ,  A]=butter(  1,0.1);  %color  the  white  noise 
n=filter(B,A,w); 


x=s+n; 

white=filter(A,B,x); 

PSDn=abs(fft(x)).A2; 

PSDw=abs(fft(white)).''2; 

P(ind,:)=PSDw(l:512); 

Pxa(ind,:)=PSDn(l:512); 

end 
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Pyave=sum(P)./numb; 

Pxave=siim(Pxa)  ./numb; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Output  resxilts 

clf 

subplot(311) 

plot((0:5 1  l)./1024,Pxave./max(Pxave)) 
ylabel(‘Averaged(abs(fft(x)).'^2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot(( 130: 1 80)/1024,Pxave( 130: 180)./max(Pxave( 130: 180))) 
axis([130/1024  180/10240  1]) 

axes(‘position’,[.7  .78  .15  .12]) 

plot((285:335)/1024,Pxave(285:335)/max(Pxave(285:335))) 
axis([285/1024  335/1024  0  1]) 

print  invfigx 

!ps3epsi  <invfigx.ps  >invfrgx.epsi 
clf 

subplot(31 1) 

plot((0:5 1  l)./l,Pyave./max(Pyave)) 
ylabel(‘Averaged(abs(fft(y)).''2)’) 
xlabel(‘Fraction  of  sampling  frequency’) 

print  invfigy 

!ps3epsi  <invfigy.ps  >mvfigy.epsi 

%%%%%%%%%%%%%%7o%%%%%%%%%%%%%%%%%%%%%%%% 
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APPENDIX  J  -  MATLAB  PROGRAM  IMPLEMENTING  THE 


EIGENVECTOR  FACTORIZATION 


%  Diagonalization  of  the  Correlation  matrix  of  the  noise 
%  ref  p50, 2.6.1 

%  do  transformation  y=E*T*x  with  E*T  of  the  noise 
%  Use  theoretical  Rn  estimated  from  PSD 


%  Author:Capt  M.A.  Cloutier 
%  Date:  2  Dec  95 
%  Name:  eigfactfin.m 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

clear 

clf 

orient  tall 
% 

%  generate  data 
numb=10; 

PSDxa=zeros(numb,5 12); 

PSDya=zeros(numb,5 12); 

for  ind2=l:numb 
1=1024; 
t=l:l; 

s=  .475*cos(2*pi*l/20)  +0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10); 

%varia=input(‘What  is  the  variance  of  the  white  noise?  ‘); 
varia=10; 

w=sqrt(varia)*randn(  1  ,length(s));%white  noise 

[B,A]=butter(l,0. l);%color  the  white  noise 
n=filter(B,A,w); 


x=s+n; 


%calculate  autocorrelation  sequence  of  the  noise 


117 


[H,P]=freqz(B,A,2048);%H  contains  the  transfer  function 


Rnt=^eal(ifft(abs([H.’  fliplr(H.  ’)])  A2));%with 
correlation  matrix 

Rnmat=toeplitz(Rnt(l:256)); 


Wiener-Khinchine  theorem  find 


[V,D]=eig(Rnmat);%  D:diagonal  matrix  of  eigenvalue, V  eigenvectors  x*V=V*D 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

index=l; 


%do  transformation  of  x  in  blocks  of  256  points 


fori=l:4 

y(index;index+255)=V’  *x(index:index+255). 
index=index+256; 
end 


%%%%%%%7o%%%%%%%%%%%%%%%%%%%%%% 
%Calculate  average  PSD  for  x  and  y 

PSDx=abs(fft(x,  1024)A2); 

PSDxa(ind2,:)=PSDx(l  :5 12); 

PSDy=abs(ffl(y, 1024))  A2; 

PSDya(ind2.:)=PSDy(l  :5 12); 

end 


PSDyave=sum(PSDya)ynumb; 

PSDxave=sum^SDxa)ynumb; 


%Ou^ut  results 


clf 

subplot(311) 

plot((0:511)/1024,PSDxave(l:512)./max(PSDxave(l:512))); 

ylabel(‘Averaged(abs(fiit(x))A2)’) 

xlabel(‘Fraction  of  sampling  frequency’) 


axes(‘position’,[.45  .78  .15  .12]) 

plot((130:180)/1024,PSDxave(130:180)./max(PSDxave(130:180))) 
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axis([130/1024  180/10240  1]) 
axes(‘position’,[.7  .78  .15  .12]) 

plot((285:335)/1024,PSDxave(285:335)/max(PSDxave(285:335))) 
axis([285/1024  335/1024  0  1]) 

print  eigfactfigx 

!ps3epsi  <eigfactfigx.ps  >eigfactfigx.epsi 
clf 

subplot(311) 

plot((0:5 1  l)/1024,PSDyave(  1 :5 12)  ./max(PSDyave(  1 :5 12))); 

ylabel(‘Averaged(abs(fft(y)).''2)’) 

xlabel( ‘Fraction  of  sampling  frequency’) 

print  eigfactfigy 

!ps3epsi  <eigfactfigy.ps  >eigfactfrgy.epsi 
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APPENDIX  K  -  MATLAB  PROGRAM  IMPLEMENTING  THE 


SINGULAR  VALUE  DECOMPOSITION 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%7o% 
%  Single  Value  Decomposition:  From  singvdbigN.m 
% 


%  Author:Capt  M.A.  Cloutier 
%Date:  15  Nov  95 
%Name:  svdtin.m 


clear 

clf 

%  generate  data 

1=1024; 

t=l:l; 


s=.475*cos(2*pi*t/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 
varia=50; 


PSDy=zeros(10,1024); 
PSDx=z«-os( 10, 1024); 
average=20; 


for  j=l:average 

w=sqrt(varia)*randn(l,l);%white  noise 

[B , A]=butter(  1,0.1);  %color  the  white  noise 

n=filter(B,A,w); 


x=s+n; 


%Defme  matrix  of  sample  vectors  of  the  noise:  Use  48  vectors  1024  pts  long 


fori=l:48 

N(i,:)=n(l:1024); 

w=sqrt(varia)*randn(l,1024);%white  noise 
n=filter(B,A,w); 
end 

%Perform  the  SVD  on  the  noise  data  matrix  N 
%  N=UEVh:  where  V  is  the  matrix  of  eigenvectors 

[U,S,V]=svd(N); 

%do  transformation  in  one  block  of  1024  points 
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index=l; 

y(index:index+1023)=V’*x(index:index+1023).’; 

PSDy(j,:)=abs(fft(y,1024))A2; 

PSDx(j,:)=abs(fft(x,1024))A2; 

end 

PSDyave=sum(PSDy(:  ,1:5 12))yaverage; 

PSDxave=sum(PSDx(:,  1 :5 12))./average; 

clf 

subplot(311) 

plot((0:5 1 1 )./1024,PSDxave./max(PSDxave)) 
ylabel(‘Averaged(abs(fft(x)).''2)’) 
xlabel(‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot((130:180)/1024,PSDxave(130:180)./max(PSDxave(130:180))) 
axis([130/1024  180/1024  0  1]) 

axes(‘position’,[.7  .78  .15  .12]) 

plot((285:335)/1024,PSDxave(285:335)/max(PSDxave(285;335))) 
axis([285/1024  335/1024  0  1]) 

print  svdfigx 

Ips3epsi  <svdfigx.ps  >svdfigx.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./lJPSDyave./max(PSDyave)) 
ylabel(‘Averaged(abs(fft(y)).''2)’) 
xlabel(‘Fraction  of  sampling  frequency’) 


print  svdfigy 

!ps3epsi  <svdfigy.ps  >svdfigy.epsi 


APPENDIX  L  -  MATLAB  PROGRAM  IMPLEMENTING  THE 
MAHALANOBIS  TRANSFORMATION 


%  Mahalanobis  Transformation  using  the  Correlation  matrix  of  the  noise 
%  ref  p247,  eq  5.63 

%  do  transformation  y=Rn''-.5*x  with  Rn  of  the  noise 


%  Author:Capt  M.A.  Cloutier 
%  Date:  5  Dec  95 
%  Name:  mahafm.m 


clear 

clf 

orient  tall 


numb=10; 

PSDx=zeros(numb,5 12); 
PSDy=zeros(numb,5 12); 


for  ind=l:numb 
%  generate  data 
1=1024; 
t=l:l; 

s=  .475*cos(2*pi*l/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10); 
varia=25; 

w=sqrt(varia)*randn(  1  ,length(s));%white  noise 
[B,A]=butter(l,0, 1);  %color  the  white  noise 
n=fllter(B,A,w); 


x=s+n; 


[H,P]=freqz(B  A,256);%H  contains  the  transfer  function 
Rnt=real(iffi(abs([H,’ fliplr(H.’)]).'^2));%with  Wiener-Khinchine  theorem  find 
correlation  matrix 

Rnmat=toeplitz(Rnt( 1 :256)); 

index=l; 

%do  transformation  in  blocks  of  256  points 
for  i=l:4 

y(index:index+255)=inv(sqrt(Rnmat))*x(index:index+255).’; 

index=index+256; 
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end 


Py=abs(fft(y,1024))A2; 

Px=abs(fft(x,1024))A2; 

PSDy(ind,:)=Py(l:512); 

PSDx(ind,:)=Px(l:512); 

end 

Pyave=sum(PSDy)./numb; 

Pxave=siini(PSDx)./niimb; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Output  results 

clf 

subplot(311) 

plot((0:5 1  l)./1024,Pxave./max(Pxave)) 

ylabel(‘Averaged(abs(fft(x)).''2)’) 

xlabel(‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot((130:180)/1024,Pxave(130:180)./max(Pxave(130:180))) 

axis([130/1024  180/10240  1]) 

axes(‘position’,[.7  .78  .15  .12]) 

plot((285:335)/1024,Pxave(285:335)/max(Pxave(285:335))) 

axis([285/1024  335/1024  0  1]) 

print  mahaflgx 

!ps3epsi  <mahafrgx.ps  >mahafigx.epsi 
clf 

subplot(311) 

plot((0;5 1 1)./1 024,Pyave./max(Pyave)) 
ylabel(‘Averaged(abs(fft(y)).'^2)’) 
xlabel(‘Fraction  of  sampling  frequency*) 
axes(‘position’,[.17  .78  .12  .12]) 
plot((25:75)/1024,Pyave(25:75)/max(Pyave(25:75))) 
axis([25/1024  75/1024  0  1]) 

print  mahafrgy 

!ps3epsi  <maha£igy.ps  >mahafrgy.epsi 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%^^^^ 
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APPENDIX  M  -  MATLAB  PROGRAM  IMPLEMENTING  THE  LDU 
AND  UDL  TRIANGULAR  FACTORIZATIONS 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Diagonalization  by  triangular  decomposition  using  LDU  and  UDL  decomp 
%  LDL:  y=inv(L)*x ;  UDL:  y=inv(Lu)*x 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  5  Dec  95 
%  Name:  LDLfin.m 

%%%%%%%7o%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  ref  p64, 2.7.1 
% 

clear 

clf 

orient  tall 
numb=10; 

PSDx=zeros(numb,5 12); 

PSDy=zeros(numb^l2); 

for  ind=l:numb 

%  generate  data 

1=1024; 

t=l:l; 

s=  .475*cos(2*pi*l/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10); 
varia=25; 

w=sqrt(varia)*randn(  1  ,length(s));%white  noise 

[BA^]=butter(l,0.1);  %color  the  white  noise 
n=filter(B,A.w); 

x=s+n; 

[H,P]=fTeqz(B,A,2047);%H  contains  the  transfer  function 

Rnt=real(ifft(abs([H.’  fliplr(H.’)]).''2));%with  Wiener-Khinchine  theorem  find 
correlation  matrix 

Rnmat=toeplitz(Rnt( 1 :25 6)); 
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[L,U,P]=lu(Rnmat);%for  LDU 

[Lu,Uu,Pu]=lu(£liplr(flipud(Rnmat)));%UDL:LDU  or  reversed  Rn 
%do  transformation  y=inv(L)x 
index=l; 
for  i=l:4 

y(index:(index+255))=inv(L)*x(index;(index+255)).’; 

yh(index:(index+255))=inv(fliplr(flipud(Lu)))*x(index:(index+255)).’; 

index=index+256; 

end 

Py=abs(fft(y,1024))A2; 

Pyh=abs(fft(yh,1024))A2; 

Px=abs(fft(x.l024)).'^2; 

PSDy(ind,:)=Py(l:512); 

PSDyh(ind,;)=Pyh(l:512); 

PSDx(ind,:)=Px(l:512); 

end 

Pyave=sum(PSDy)./numb; 

Pyaveh=sum(PSDyh)./numb; 

Pxave=sum(PSDx)./numb; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Ou^ut  results 

clf 

subplot(311) 

plot((0:5 1 1 )./1024JPxave./max(Pxave)) 
ylabel(‘Averaged(abs(fft(x)).''2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot((130: 180)/1024,Pxave(130: 180)./max(Pxave(130: 180))) 

axis([130/1024  180/10240  1]) 
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axes(‘position’,[.7  .78  .15  .12]) 
plot((285:335)/1024,Pxave(285:335)/max(Pxave(285:335))) 
axis([285/1024  335/10240  1]) 

print  Idlfigx 

!ps3epsi  <ldlfigx.ps  >ldlfigx.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./1024JPyave./max(Pyave)) 
ylabel(‘Averaged(abs(fft(y)).'^2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  Idlfigy 

!ps3epsi  <ldlfigy.ps  >ldlfigy.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./1024,Pyaveh./max(Pyaveh)) 
ylabel(‘Averaged(abs(fft^)).A2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  Idlhfrgy 

!ps3epsi  <l(flhfigy.ps  >ldlhfigy.epsi 
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APPENDIX  N  -  MATLAB  PROGRAM  IMPLEMENTING  THE 


CHOLESKY  FACTORIZATION 


%%%%%%%%%%%%%%%%%%%%7o%%%%%%%%%%%%%%%% 

%  Diagonalization  by  triangular  decomposition;  Cholesky 
%  y=inv(L)*x  (upper)  and  y=inv(fliplr(flipud(Ru’)))*x  (lower) 

% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  5  Dec  95 
%  Name:  cholfin.m 

%%%%%%%%%%%%%%%%7o%%%%%%%%%%%%%%%%%%%% 

% 

clear 

clf 

orient  tall 
numb=10; 

PSDx=z«ros(munb,5 12); 

PSDy=zeros(numb,5 12); 

for  ind=l:numb 
%  generate  data 
1=1024; 
t=l:l; 

s=  .475*cos(2*pi*t/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 
varia=25; 

w=sqrt(varia)*randn(l,length(s));%white  noise 

[B ,  A]=butter(  1 ,0- 1  )i  %color  the  white  noise 
n=filter(B,A,w); 


x=s+n; 

[H,P]=freqz(B  A^,2047);%H  contains  the  transfer  function 

Rnt=real(ifft(abs([H.’ fliplr(H.’)])A2));%with  Wiener-Khinchine  theorem  find 
correlation  matrix 

Rnmat=toeplitz(Rnt( 1 :256)); 

R=chol(Rnmat);  %  for  lower-upper 
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Ru=chol(fliplr(flipud(Rnmat)));%  for  upper  lower 
%do  transformation  y=inv(R’)x 


index=l; 

%The  matlab  chol  produces  an  upper  triangular 

for  i=l:4 

y(index:(index+255))=inv(R’)*x(index:(index+255)).’; 
yh(index:(index+255))=inv(fIiplr(flipud(Ru’)))*x(index:(index+255)).’; 
index=index+25 6; 
end 


Py=abs(fft(y,1024))A2; 
Pyh=abs(fft(yh, 1024)).'^2; 
Px=abs(fft(x,1024))A2; 

PSDy(ind,:)=Py(l:512); 

PSDyh(ind,:)=Pyh(l:512); 

PSDx(ind,:)=Px(l:512); 

end 


Pyave=sum(PSDy)./numb; 

I^aveh=sum(PSDyh)./numb; 

Pxave=sum(PSDx)./numb; 


%  Output  results 
clf 

subplot(31 1) 

plot((0:5 1 1 )  ./1024,Pxave./max(Pxave)) 
ylabel(  ‘  Averaged(abs(fft(x))  A2)’ ) 
xlabel(‘Fraction  of  sampling  frequency’) 


axes(‘position’,[.45  .78  .15  .12]) 
plot((130:180)/1024,Pxave(130:180)./max(Pxave(130:180))) 
axis([ 130/1024  180/10240  1]) 


axes(‘position’,[.7  .78  .15  .12]) 
plot((285:335)/1024,Pxave(285:335)/max(Pxave(285:335))) 
axis([285/1024  335/1024  0  1]) 
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print  chotQgx 

!ps3epsi  <cholAgx.ps  >cholfigx.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./1024,Pyave./niax(Pyave)) 
ylabel(  ‘  Averaged(abs  (ff t(y ))  .'^2)  ’ ) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  cholfigy 

!ps3epsi  <cholfigy.ps  >cholfigy.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./1024,Pyaveh./max(Pyaveh)) 
ylabel(‘Averaged(abs(fft(y))A2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  cholhfigy 

!ps3epsi  <cholhfigy.ps  >cholhfigy.epsi 


APPENDIX  O  -  MATLAB  PROGRAM  IMPLEMENTING  THE  QR 

FACTORIZATION 


%  Diagonalization  by  triangular  decomposition 
% 

%  QR;  y=dnv(L)*x 

%  QR  uses  a  data  matrixrThe  QR  factorization  provides  the  factors  in  the 
%  triangular  decomposition  of  the  correlation  matrix. 

%  X=QR:  data  matrix  X  is  expressed  as  the  product  of  a  rectangular  matrix 
%  whose  columns  are  orthonormal  and  a  square  upper  triangular  matrix, 

% 

%  This  version  uses  a  256X256  noise  data  matrix 


%  Author:Capt  M.A.  Cloutier 
%  Date:  26  Oct  95 
%  Name:  QRbigfin.m 


% 


clear 

clf 

orient  tall 


numb=10; 

PSDx=zeros(numb,5 12); 
PSDy=zeros(numb,5 12); 

for  ind=l:numb 


%  generate  data 

1=1024; 

t=l:l; 

s=  .475*cos(2*pi*l/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 
varia=50; 


w=sqrt(varia)*randn(l,256*256);%white  noise 

[B , A]=butter(  1,0.1);  %color  the  white  noise 
n=filter(B,A,w); 

x=s+n(l:1024); 
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%Defme  matrix  of  sample  vectors  of  the  noise:  Use  256  vectors  256  pts  long 

N=zeros(256,256); 

index=l; 

fori=l:256 

N(i,:)=n(index:(index+255)); 

index=index+256; 

end 

[Q.R]=qr(N); 

Dhalf=(l/sqrt(256))*diag(diag(R)); 

L=(l/sqrt(256))*R’*inv(Dhalf); 

%do  transformation  y=dnv(L)x 

%do  transformation  in  blocks  of  256  points 

index  =1; 

fori=l:4 

y(index:index+255)=inv(L)*x(index:index+255).’; 

index=dndex+256; 

end 

Py=abs(fft(y,1024)).''2; 

Px=abs(fft(x,1024))A2; 

PSDy(ind,:)=Py(l:512); 

PSDx(ind,:)=Px(l:512); 

end 

Pyave=sum(PSDy)./numb; 

Pxave=sum(PSDx)./numb; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Output  results 

clf 

subplot(311) 

plot((0:5 11)  ./1024,Pxave./max(Pxave)) 
ylabel(‘Averaged(abs(fft(x)).'^2)’) 
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xlabel( ‘Fraction  of  sampling  frequency’) 


axes(‘position’,[.45  .78  .15  .12]) 

plot((130: 180)/1024,Pxave(130: 180)./max(Pxave(130: 180))) 
axis([ 130/1024  180/10240  1]) 

axes(‘position’,[.7  .78  .15  .12]) 
plot((285:335)/1024,Pxave(285:335)/max(Pxave(285:335))) 
axis([285/1024  335/1024  0  1]) 

print  QR2figx 

!ps3epsi  <QR2figx.ps  >QR2figx.epsi 
clf 

subplot(311) 

plot((0:5 1  l)./1024,Pyave./max(Pyave)) 
ylabel(‘Averaged(abs(fft(y)).^2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  QR2figy 

!ps3epsi  <QR2figy.ps  >QR2figy.epsi 
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APPENDIX  P  -  MATLAB  PROGRAM  IMPLEMENTING  THE 


DIFFERENTIAL  WHITENING  FILTER 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Noise  Whitener  based  on  differentiation 
% 

%  Author:Capt  M.A.  Cloutier 
%  Date:  29  Nov  95 
%  Name:  differefin.m 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear 

clf 

orient  tall 
numb=10; 

PSDx=zeros(numb,5 12); 

PSDy=zeros(numb,5 12); 

PSDy2=zeros(numb,5 12); 

PSDy3=zeros(numb,5 12); 


%%%%%%%%%%%%%%%%%%%%%%% 
for  ind=l:numb 

1=1024; 

t=l:l; 

s=  .475*cos(2*pi*t/20)  +  0.2*cos(2*pi*t*3/20)  +  0.07*cos(2*pi*t*3/10  +  pi/3); 
varia=25; 

w=sqrt(varia)*randn(  1  ,length(s)) ;  %white  noise 

[B,A]=butter(l,0. 1);  %color  the  white  noise 
n=filto'(B,A,w); 

x=s+n; 


%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  First  order  differentiation 

y=diff(x); 
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Py=abs(fft(y,1024)).''2; 

Px=abs(fft(x,1024)).'^2; 


PSDy(ind,:)=Py(l:512); 

PSDx(ind,:)=Px(l:512); 

7o%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Second  order  differentiation 

y2=diff(x.2); 

Py2=abs(fft(y2,1024))A2; 

PSDy2(ind,:)=Py2(l:512); 


%%%%%%%%%%%%%%%%7o%%%%%%%%%% 
%  Third  order  differentiation 

y3=diff(x,3): 

Py3=abs(fft(y3, 1024))  A2; 

PSDy3(ind,:)=Py3(l:512); 

end 

Pyave=sum(PSDy)./nunib; 

Pxave=sum(PSDx)./numb; 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Output  results 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  First  order 

clf 

subplot(31 1) 

plot((0:5 1  l)./1024,Pxave./niax(Pxave)) 
ylabel(‘Averaged(abs(fft(x))  A2)’) 
xlabel(‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.45  .78  .15  .12]) 

plot((130: 180)71024, Pxave(130: 180)./max(Pxave(130: 180))) 
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axis([ 130/1024  180/10240  1]) 

axes(‘position’,[.7  .78  .15  .12]) 
plot((285:335)/1024,Pxave(285:335)/max(Pxave(285;335))) 
axis([285/1024  335/1024  0  1]) 

print  difffigx 

!ps3epsi  <difffigx.ps  >difffigx.epsi 
clf 

subplot(311) 

plot((0:5 1 1)./1024 J*yave./max(Pyave)) 
ylabel(‘Averaged(abs(fft(y)).''2)*) 
xlabel( ‘Fraction  of  sampling  frequency’) 

print  difffigy 

!ps3epsi  <difffigy.ps  >difffigy.epsi 


%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Second  order 

Pyave2=sum(PSDy2)./numb; 

clf 

subplot(311) 

plot((0:5 1  l)./1024J’yave2./max(Pyave2)) 
ylabel(‘Averaged(abs(fft(y)).''2)’) 
xlabel( ‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.17  .78  .12  .12]) 
plot((25:75)/1024,Pyave2(25:75)/max(Pyave2(25:75))) 
axis([25/1024  75/1024  0  1]) 

print  diff2figy 

!ps3epsi  <diff2figy.ps  >diff2figy.epsi 

%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Third  order 

Pyave3=smn(PSDy3)./numb; 

clf 
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subplot(311) 

plot((0:5 1  l)./1024,Pyave3./max(Pyave3)) 
ylabel(  ‘  Averaged(abs(fft(y))  .'^2)  ’ ) 
xlabel( ‘Fraction  of  sampling  frequency’) 

axes(‘position’,[.17  .78  .12  .12]) 
plot((25:75)/1024,Pyave3(25:75)/max(Pyave3(25:75))) 
axis([25/1024  75/102401]) 

print  diffSfigy 

!ps3epsi  <diff3figy.ps  xiifBfigy.epsi 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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APPENDIX  Q  -  MATLAB  PROGRAM  IMPLEMENTING  THE 
WHITENING  PREDICTION  ERROR  FILTER 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Whitening  property  of  prediction  error  filter 


%  Author:  Capt  M. A.  Cloutier 
%Date:  16  Nov  95 
%Name:  adaptfm.m 


clear 

figure(l) 

clf 

(Hienttall 

1=18432;  %  length  of  the  data 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%Generate  the  data 
% 

t=l:l; 

s=  .475*cos(2*pi*t/8)  +0.15*cos(2*pi*t*2/8)  +  0.1*cos(2*pi*t*3/8  +  pi/3);% 
varia=100;  %  variance  of  the  AWGN  to  be  coltx’ed 


w=sqrt(varia)*randn(  1  ,length(s));%white  noise 

[B,A]=butter(l,0.1);%color  the  white  noise  Butter  1st  order 
n=filter(B,A,w); 


X  =  s+n;  %  X  is  the  signal  +  noise 

a=input(‘What  filter  size  is  required  ?’); 
mu=0.00(X)5;  %size  of  the  learning  parameter 


w  =  ones(l,a);  %  Initialize  filter  weights  to  1 


X  =  [zeros(l,a-l)^];  %pad  input  to  start  data 
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nest  =  zeros(  1  ,l);%initialize  the  predicted  error  vector 

skest  =  zeros(  1  ,l);%initialize  the  estimated  signal  vector 

%Estimate  the  correlation  matrix  of  the  noise  based  on  its  PSD 

[H,P]=freqz(B,A,2047);%H  contains  the  transfer  function 

Rnt=abs(ifft(abs(H) A2));%with  Wiener-Khinchine  theorem  find  Rn 

Rs=xcorr(s(l;1024))./1024;%calculate  autocorrelation  of  signal 


%%%%%%%%%%%%%%7o%%%%%%%%%%%%% 

%Plot  both  autocorrelation  functions  to  identify  maxs  and  zero  crossings 
clf 

subplot(311) 
plot(0:20,  Rnt(l:21)) 
grid 

ylabel(‘Amplitude’) 
xlabel(*Delay’) 
print  adapfigRn 

!ps3epsi  <adapfigRn.ps  >adapfigRn.epsi 
clf 

subplot(311) 

plot((0:20),Rs(1024: 1044)) 
grid 

xlabel(‘Delay’) 
ylabel(‘Amplitude’) 
print  adapfigRs 

!ps3epsi  <a^pfigRs.ps  >adapfigRs.epsi 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Rim  adaptive  LMS  algorithm  on  all  data 

delay=input(‘What  delay?  ‘);%choose  delay  based  on  correlation  functions 
for  i=l:l-delay 

skest(i)  =  w*flipud(x(i;i+(a-l))’);%  compute  signal  est 
nest(i)  =  x(i+a+delay-l)-skest(i);%  compute  predicted  errcH- 
w  =  w+mu.*nest(i).*fliplr(x(i'(i+(a-l)))); 
end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Compute  averaged  peiiodograms  for  the  input  signal,  predicted  error 
%  and  estimated  signal. 

PSDl=zeros(10,1024); 

PSD2=zeros( 10,1024); 

PSD3=zeros( 10,1024); 

index=8193; 

fori=l:10 

PSDl(i,:)=abs(fft(x(index;index+1023),1024)).A2; 
PSD2(i,:)=abs(fft(nest(index:mdex+1023),1024)).''2; 
PSD3(i,:)=abs(fft(skest(index:index+1023),1024)).''2; 
index=index+ 1024; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Output  average  PSDs 

clf 

subplot(311) 

PSDx=sum(PSDl)./10; 

PSDx(  1 :5 12)=PSDx(l  :5 12)./max(PSDx(l  :5 12)); 
plot((0:51  l)./1024,PSDx(l;5 12)) 
ylabel(‘Averaged(abs(fft(x)).'^2)’) 
xlabel( Traction  of  sampling  frequency’) 
axis([0 .5  0  1]) 

axes(‘position’,[.59  .78  .3  .12]) 
plot((241;400)/1024,PSDx(241:400)./max(PSDx(241:400))) 
axis([241/1024  400/10240  1]) 

print  adapfigx 

!ps3epsi  cadapfrgx.ps  >adapfigx.epsi 
clf 

subplot(311) 

PSDer=sum(PSD2)./10; 

plot((0:5 1  l)./1024,PSDer(l  :5 12)./inax(PSDer(l  :5 12))) 

ylabel(‘Averaged(abs(fft(x)).''2)’) 

xlabel( ‘Fraction  of  sampling  frequency’) 

axis([0 .5  0  1]) 

print  adapfigerr 

!ps3epsi  <adapfigerr.p8  >adapfrgerr.epsi 


143 


subplot(31 1) 

PSDsig=sum(PSD3)./10; 

plot((0:5 1  l)./1024,PSDsig(  1 :5 12)./max(PSDsig(  1 :5 12))) 
ylabel(  ‘  Averaged(abs(fft(x)).'^2)  ’ ) 
xlabel(‘Fraction  of  sampling  frequency’) 
axis([0 .5  0  1]) 

axes(‘position’,[.59  .78  .3  .12]) 

plot((241:400)/1024,PSDsig(241:400)./max(PSDsig(241:400))) 
axis([241/1024  400/1024  0  1]) 


print  adapfrgsh 

!ps3epsi  <adapfigsh.ps  >adapfrgsh.epsi 
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